A mathematical modelling approach to uncover factors influencing the spread of Campylobacter in a flock of chickens

Despite continued efforts into improving biosecurity protocol, Campylobacter continues to be detected in the majority of commercial chicken flocks across Europe. Using an extensive data set of Campylobacter prevalence within a chicken breeder flock for over a year, multiple Bayesian models are presented to explore the dynamics of the spread of Campylobacter in response to seasonal variation, species-specificity, bird health and total infection prevalence. It was found that birds within the flock varied greatly in their response to bacterial challenge, and that this phenomena had a large impact in the overall prevalence of different species of Campylobacter. Campylobacter jejuni appeared more frequently in the summer, while Campylobacter coli persisted for a longer duration, amplified by the most susceptible birds in the flock. Our study suggests that strains of Campylobacter that appear most frequently likely possess no demographic advantage, but are instead amplified due to the health of the birds that ingest it.

the flock. This is shown below in Figure 1, with all positive samples classified by species of Campylobacter.  Within each species, multiple STs are recorded. In figures 2 and 3 below we plot the five week moving averages of total positive 87 samples for each species. Beneath each point we plot a histogram showing how this average is split between the competing STs. 88 We notice from figures 2 and 3 that there are more significant STs of C. jejuni than C. coli, despite both species existing at 89 roughly equal levels. We also see that C. jejuni appears to peak in the summer, around the August period, coinciding with a dip 90 in the population of C. coli STs. Within each species we can observe that different ST populations grow and shrink across 91 the study period. For example, within figure 2 we see that the summer peak is dominated by the prevalence of ST 51 and 53, 92 however by November/December, this population shrinks, and instead ST 607 rapidly increases in population.  STs that appear less than ten times throughout the entire experiment are amalgamated into a group "Low Count".

5/33
through querying the probability of chickens transitioning from different infection states using a series of Bayesian models 98 presented below. 99 Model Development 100 In this section we discuss the general methodology behind all of our models. A general step-by-step process to model formulation is also presented in Box 1. Each model begins by classifying each of the datapoints into certain state labels. For example, at the simplest level each reading can be classified as either "State 1: Uninfected" or "State 2: Infected". Other models may use more states to further distinguish infections by species or ST. After doing this, we are able to convey this classification data in the form of a matrix S[c,t] where c ∈ {1, 2, ..., 200} is the index denoting which chicken is considered, and t ∈ {1, 2, ..., 51} is the index denoting which week is considered. Therefore each element of S will be a number conveying the state classification of that particular data point. For example, S[3, 7] = 1, would indicate that on week 7, chicken number 3 was classified as state 1; uninfected. Because only 75 of the 200 chickens were tested at random each week, many of these matrix elements are undefined, and as such are marked as 'NA'.
Once the matrix is defined, each model uses a Bayesian process to find the transition probabilities between these states. In short, π i, j is the probability that a chicken moves from state i to state j across a week. The exact choice of how to formulate the expressions is where our models vary, as different formulations are able to investigate different relationships governing these transition probabilities. For example, at the simplest level, we could define where we seek to find the values α 1 ∈ [0, 1] and α 2 ∈ [0, 1] that best fit the data S. Note that we have bounded π i, j between 0 and 1, as each value represents a probability. Likewise each row of π must sum to 1, as these probabilities cover all transition possibilities. In the example of equations (1) above, when starting from state 1, one can transition to state 2 (π 1,2 ), or remain in state 1 (π 1,1 ), hence π 1,1 + π 1,1 = 1. Different models below will use more complex definitions for π to explore the impact of time, density dependence, and chicken health on transitions between different states. A Bayesian statistical model provides a way to iteratively deduce parameters of interest in regards to given data. The

6/33
process is derived from Bayes' theorem: where θ is the parameter/s we wish to discover, and D is the data provided. In short, equation (2) reads that when starting from 101 an initial, prior, belief in what values θ may take (P(θ )), one may obtain an updated, posterior, probability distribution for 102 these possible values given some provided data (P(θ |D)). A more thorough introduction to Bayesian modelling is provided in 103 Appendix 1. In our case, the parameters we seek, θ , are the ones used in our definition of π, such as α 1 and α 2 in the example 104 above. The data, D, we use is the matrix S. Choose how data should be classified, and construct matrix S containing all state classifications for each data point.

Decide formulation of transition matrix.
Choose how model will define transition probabilities and dependencies.

Run Bayesian model.
Define prior probability distributions for model parameters. Program and run Bayesian model using JAGS, to acquire a posterior probability distribution for all model parameters defined in step 2.

Assess convergence.
Investigate model output to assure posterior distribution is well-constructed and has converged.

Present results.
Plot the transition probabilities, π i, j , and interpret the results.

107
Below we present a series of case studies presenting our different models and their results. All models were run using JAGS 31 108 from within R using the run.jags package 32 . Our first model investigates how time affects the transition probabilities between states. Following the process outlined in Box 1, we choose to initially classify our data as one of two states: "state 1: uninfected" and "state 2: infected".
To assess how the transition probabilities vary through time we must ensure that we define our transition probabilities such that they depend on time. One way would be to adapt equations (1) above such that π 1,1 was a function of α + βt.
However, this would impose structure upon the transition probabilities, enforcing them to change linearly with time. Ideally

7/33
a model formulation should allow as much freedom as possible to fit to the data. As such, we shall instead construct π as a three-dimensional array. In essence this means that each time period can be described by its own transition matrix. Formally we write this as, for t ∈ {1, 2, ..., 51}. Here ilogit() is the inverse logit function defined by ilogit(x) = e x 1+e x . This function is bounded between 0 The results for this model are presented below in figure 4. The median values of the transition probabilities for (4A) 137 π 1,1,t , (4B) π 1,2,t , (4C) π 2,1,t , (4D) π 2,2,t are plotted, and a linear regression is fit to these outputs using the lm function in R.

138
The probability of transitioning from state 1 (plots 4A and 4B) was not significantly correlated against time (t-test, p = 0.135), 139 however transitions from state 2 (plots 4C and 4D) against time were statistically significant (t-test, p < 0.01).
140 Figure 4. Transition probabilities between two states, 'uninfected' and 'infected'. Plots show (A) π 1,1,t , (B) π 1,2,t , (C) π 2,1,t and (D) π 2,2,t against time. Each point is the calculated transition probability for that time point. Also plotted is a linear regression against these points in blue, with a shaded region depicting the 95% confidence interval of the regression. (C) and (D) are significant (p < 0.01).
These findings suggest that infected chickens were more likely to remain infected, and less likely to clear infection, as time 141 progressed.
The model was run with two chains and an initial burn-in period of 5,000 iterations. Posterior distributions were built from a  each state are very similar, however one can note that a chicken is more likely to be infected with C. coli when transitioning 147 from a state of already being infected with C. coli. We also see that a chicken infected with C. coli is less likely to transition to 148 being uninfected than a chicken infected with C. jejuni.

149
Model 3: Time and species dependence 150 We now combine the previous two models together, to investigate how the transitions between species alter across time. We once again therefore classify our data into three categories, as per the previous model.
We will be constructing a three-dimensional array once again for our transition probabilities, with each time period be-10/33 ing described by a separate 3 × 3 transition matrix. To ensure each row of these matrices sums to 1, we start by framing the transition probabilities as an unbounded array p, before scaling these into our final array π. p is defined as The exponential function here assures that, like in our initial model, our α parameters will describe the average transition value across time, with the C parameters describing a small perturbation away from this mean. C values only need to be implemented on two probabilities in each row, as we will next scale these so that each row sums to 1, meaning that two free correction terms are sufficient to describe the distribution of the row. Our scaling is then performed like so, π 2,2,t = p 2,2,t p 2,1,t + p 2,2,t + p 2,3,t , π 2,3,t = p 2,3,t p 2,1,t + p 2,2,t + p 2,3,t , π 3,1,t = p 3,1,t p 3,1,t + p 3,2,t + p 3,3,t , π 3,2,t = p 3,2,t p 3,1,t + p 3,2,t + p 3,3,t , π 3,3,t = p 3,3,t p 3,1,t + p 3,2,t + p 3,3,t .
We choose priors of N(0, 1000) for all our α values (normal distributions with mean 0 and standard deviation 1000). Like the 151 first model, we shall construct a hierarchical dependency such that our C i [t] are all drawn from a normal distribution for each t.

152
Motivated by the correlation observed in the first model, we actually set these six C i terms to all be drawn from a 6-variable 153 multivariate normal distribution, with mean (0, 0, 0, 0, 0, 0) and a covariance matrix as our parameter to be defined. JAGS 154 requires the input of a precision matrix (the inverse of the covariance matrix) for its formulation of the multivariate normal 155 distribution, so we set a prior distribution on the precision matrix of Wishart(I 6 , 6), where I 6 is the 6 × 6 identity matrix.
the data in figure 6, we also tested for statistical significance against a quadratic regression. A quadratic fit would be a strong 163 argument for the existence of seasonal variation, by capturing a difference in the middle of the time series as the time axis 164 moves to summer, before returning to winter. Recall again that this time period plotted is in weeks from February to February.

165
Only one transition probability was found to be statistically significant however, the transition from infection with C. jejuni to C. 166 coli, π 2,3,t (t-test, p < 0.05). This quadratic regression is presented in figure 7 below. This would correlate with the behaviour 167 observed in figures 2 and 3, whereby C. jejuni appears to be most prevalent in the summer, and C. coli most prevalent in the 168 winter. Each point is the calculated transition probability for that time point. Also plotted is a quadratic regression against these points in blue, with a shaded region depicting the 95% confidence interval. The transition probability was found to be statistically significant for correlation against time (t-test, p < 0.05).
The model was run with 2 chains for a burn-in period of 5,000 iterations before building posteriors from a final sample of

173
The most notable difference is seen in the perseverance of C. coli STs compared to C. jejuni STs. Comparing columns 3 and 5 'S1: uninfected' or 'S2: infected'. We then, like model 1, consider some average transition probability that each chicken is close to, and then consider some small "correction term" unique to each chicken, which may make them more or less likely to transition to a certain state. Formally, we write, π 1,1,t = 1 − π 1,2,t , π 2,1,t = ilogit(α 2 +C 2 [c]), π 2,2,t = 1 − π 2,1,t , The model was run with two chains for an initial burn-in period of 20,000 iterations, before posteriors were then con-185 structed from a sample of 50,000 iterations. Convergence was well-achieved, with all chains well-mixed and all parameters 186 sampled with a high ESS and MCSE < 0.01. The mpsrf was unable to be calculated due to the high number of stochastic nodes, 187 however there were no signs to suggest invalid convergence.

189
Upon calculating our transition probabilities for each bird, we plot the values for π 1,2 against the value of π 2,1 for each 190 bird and investigate the correlation. Figure 9 shows these results overlaid with a contour of the associated multivariate normal 191 distribution, indicating the probability density of the transition probabilities for the flock.
192 Figure 9. Transition probabilities for each bird in the flock from a state of infected to uninfected (y-axis) against the transition probability from uninfected to infected (x-axis). Contours show the fit of a multivariate normal distribution to the output. 193

15/33
The strong linear relation observed reveals the presence of distinct sub-groups within the flock of birds who are infected often, 194 and those who are infected very rarely.

195
Model 6: chicken and species dependence 196 We now alter the previous model to consider the differences in transition between species of Campylobacter across all birds.

197
As such, the data is instead classified into the three states: 'state 1: uninfected', 'state 2: infected with C. jejuni' and 'state 3: infected with C. coli. This model is formulated the same way as in model 3 above. The transition probabilities follow 199 the same structure as equations (4)  for a specific chicken. Plots 10A to 10C use π 1,1,c , the transition from uninfected to uninfected as the y-axis. This acts as a 212 rough metric for "bird resilience to infection", as the more resistant birds are more likely to continue being uninfected. As such 213 plots 10A to 10C depict how transitions related to each species vary according to host bird susceptibility. Plot 10D uses π 3,3,c , 214 the transition from C. coli to C. coli as the y-axis, to compare how the perseverance of C. coli affects the infection ability of C.

215
jejuni. Linear regressions are fit to all plots in figure 10, and all were found to be statistically significant (t-test, p < 0.0001).

216
It is interesting to note that the gradient of the lines in each plot are distinctly different from one another, highlighting how 217 each species responds differently to variations in host bird health. The transition probabilities of J-to-J and U-to-J against the transition probability of C-to-C.
The model formulation is then as follows, where N t is the number of birds that data is available for at time t. Here, as with previous models, α i represents some mean transition probability that all birds are clustered around, and C i [c] represents the slight correction for each bird c. Recall that the matrix S is populated by elements '1' denoting uninfected and '2' denoting infected. Therefore the expression S[i,t] − 1 for every i and t shifts this to instead be captured as '0' signifying uninfected, and '1' signifying infected. Therefore, the expression ∑ 51 i=1 S[i,t] − 1 will be a tally of exactly how many birds are recorded as being infected at time t. Therefore, the expression conveys the exact proportion of how many birds are currently infected. Note the use of N t as, for most weeks 75

17/33
birds are recorded for every t, however, as can be seen figure 1, occasionally a few more or less were recorded each week. Note however, that during the Bayesian modelling process, values for each element of S will be imputed in the process, meaning that we can choose to measure our density dependence using either just the provided data, or also the imputed data. There are merits to both approaches, and so results are included for both below. Here β i then represents our parameters signifying the strength of the density dependent effect.
The model was initialised with prior distributions of N(0,1000) for all α i and β i parameters. The chicken corrections terms C i [c] were, like above, drawn from a multivariate normal distribution of mean (0,0) whose precision matrix we seek. The precision matrix was initialised with a prior distribution of Wishart(I 2 , 2) where I 2 is the 2 × 2 identity matrix. The model was run with two chains for an initial burn-in period of 6,000 iterations and then posterior distributions built from a sample of 25,000 iterations. This was done twice with two variations of the model. One where density dependence is calculated from provided data, and one with the addition of imputed data. The posterior distributions of our model parameters were used to simulate the transition probabilities for each flock across a full range of total flock prevalences. i.e. using the median values for α i , β i , C i and the precision matrix, we are able to build the functions that the more resilient birds truly are less likely to become infected, as opposed to just never becoming exposed to particularly 226 virulent STs. Of interest here is that the probability of clearing infection (transitioning to uninfected) is affected far more by 227 flock prevalence proportion than the probability of becoming infected. Model 2 begins to investigate the differences in transition probabilities between the two species of Campylobacter ob-266 served. A three-state system of 'uninfected', 'infected with C. jejuni' and 'infected with C. coli' was considered, and the This brief section aims to convey the basic principles of Bayesian statistics, and familiarise the reader with the terminology that 504 is be used throughout the manuscript. For an in-depth explanation, I recommend the text by Kruschke (2014) 34 .

506
Bayesian statistics is derived wholly from the relationship defined by Bayes' theorem, If we consider θ as some statistical parameter we wish to infer, and D as some data informing the parameter, then equation

508
(1) expresses that the probability distribution for our value of θ , given our dataset (P(θ |D)), is proportional to the likelihood of 509 such data (P(D|θ )) multiplied by the probability distribution of θ free of any data (P(θ )). if we wished to calculate the probability that a flipped coin will land heads up, we may have a prior belief that the coin is fair.

514
However, upon observing a data set of 5 coin flips, all of which produced heads, we may update our posterior belief to reflect 515 that the coin may be biased.

517
The analytical difficulty in this calculation lies in computing P(D) = P(D|θ )P(θ )dθ , which is often near impossible 518 for realistically complex models. Fortunately modern computing power enables us to efficiently estimate our posterior distribu-519 tions through algorithms such as Gibbs sampling and other Metropolis-Hastings schemes.

521
Hierarchical systems represent multi-variable models where some parameters depend on other parameters. Returning to 522 the example of a coin flip, say the probability of heads (θ ) is dependent on the factory in which the coin was minted. The 523 probability that a coin was from a certain factory (ω) will then inform our value of (θ ). Expressed mathematically, equation ( θ , via our model formulation. As such, when provided with data on coin flips from multiple coins from different factories, we 527 obtain a posterior probability distribution of which factory a coin has come from, and the resulting probability of a coin flip 528 resulting in heads. This structure of conditional independence means that data relating specifically to one parameter can still 529 help inform the posterior of all other dependent variables, a key advantage of Bayesian inference. identical broilers. Figure A1 shows that, as expected, all strains perform equally well and are equally represented in the 536 amount being shed into the environment. Figure A2 instead shows the same model of five demographically identical strains of 537 Campylobacter within a flock of 400 birds whose strength of immune response is drawn from a normal distributed centred 538 around the value used for Figure 8. Figure A2E shows how five demographically equal strains can be sustained at broadly 539 different levels across the flock due only to variation in bird immune response. This is caused by random chance, in that 540 whichever strain is initially picked up by a super-shedder, such as the one shown in Figure A2D then sheds large amounts of 541 that strain of Campylobacter into the environment, increasing the likelihood of then infecting other birds in the flock. This 542 result greatly implies that the results shown in the data, whereby some STs seem to persist at higher levels than others in the 543 flock, is likely due to the variation in bird transition probabilities, as opposed to phenotypic differences between STs.