A multi-scale epidemic model of Salmonella infection with heterogeneous shedding

Salmonella strains colonize the digestive tract of farm livestock, such as chickens or pigs, without affecting them, and potentially infect food products, representing a threat for human health ranging from food poisoning to typhoid fever. It has been shown that the ability to excrete the pathogen in the environment and contaminate other animals is variable. This heterogeneity in pathogen carriage and shedding results from interactions between the host’s immune response, the pathogen and the commensal intestinal microbiota. In this paper we propose a novel generic multiscale modeling framework of heterogeneous pathogen transmission in an animal population. At the intra-host level, the model describes the interaction between the commensal microbiota, the pathogen and the inflammatory response. Random fluctuations in the ecological dynamics of the individual microbiota and transmission at between-host scale are added to obtain a drift-diffusion PDE model of the pathogen distribution at the population level. The model is further extended to represent transmission between several populations. The asymptotic behavior as well as the impact of control strategies including cleaning and antimicrobial administration are investigated through numerical simulation.


Introduction
Salmonella infection is the most common vector of collective food poisoning in the developed world. As such, deciphering the mechanisms of infection in humans and animals is a fundamental step towards the design of efficient epidemiological policies, in order to reduce the burden on agrifood industry and healthcare systems resulting from Salmonella zoonoses. Salmonella is a bacterial genus composed of various pathogenic strains, that colonize and infect the digestive tract of farm livestock, such as chickens or pigs, representing a threat for human health ranging from food poisoning to typhoid fever. In this type of animal epidemic, it has been shown that the ability to excrete the pathogen in the environment (e.g. in water, food, or feces) and thus to contaminate other individuals in the farm varies from one individual to another. Some of them, called super-shedders, are permanent carriers and highly excreting the pathogen without being affected, thus causing most of the epidemic spread. This super-shedder phenotype, still poorly defined, is most probably the result of interactions between the host's immune response, the pathogen and the commensal intestinal flora. Recent studies [1] have identified key interactions between Salmonella Typhimurium and its host during infection. In the host's gut, pathogenic virulence factors promote the inflammation of the epithelium. The inflammatory process modifies the nutritional environment of the gut lumen which disturbs the ecological equilibrium of the microbiota and creates new niches that are targeted by the pathogen. The induced inflammation provides to the pathogen a competitive advantage, which can be sufficient to promote its proliferation and its transmission to another host. The aim of this paper is to propose a generic multiscale modeling framework of heterogeneous pathogen transmission in a livestock, accounting for the interaction dynamics between the commensal microbiota, the pathogen and the inflammatory response at the intra-host level and transmission at between-host scale in a single animal population. This model is further extended at the metapopulation level, to model transmission between several populations. Note that our model distinguishes from the ones proposed by Zongo et al. [4][5][6] in the last decade which are based on individual-based model and age-structured models for salmonella infection.
The outline of the paper is the following. The single population model is described in section 1. The different scales and their interconnection are detailed. We start from an ODE equation modeling the intra-host evolution of the pathogen load with respect to time, including the host response. We next introduce a stochastic perturbation to this dynamics to account for biological variability. This SDE is used to derive a drift-diffusion PDE describing the evolution of a population density with respect to time and to pathogen load. The well-posedness of the model and asymptotic convergence towards a steady state is analyzed and an extension to the case with transmission of salmonella within the population through a reservoir is also proposed, leading to the coupling between a drift-diffusion PDE for the population and an ODE for the reservoir variable. Section 2 is devoted to some simulations of the two previous models, with or without a reservoir, and we observe that in both cases, the solution always converges towards a steady state. In section 3, some simulations of epidemic control strategies based on cleaning or drug treatment, are performed and compared. Cleaning is modeled through the addition of a term in the ODE for the reservoir variable, whereas drug treatment is described by the addition of a drift term in the population PDE equation. Then, in section 4, a compartment model taking into account some exchanges between populations is introduced and studied numerically. The model becomes in this section a large system that couples the previous PDE-ODE model, defined in each compartment, through transfer flux terms. Finally, some conclusive remarks and perspectives are given.

Mathematical models for a population structured by pathogen load
In this section, we will first present an ODE model which describes the evolution of pathogens in the gut for an individual. From this model, we will deduce a PDE for a population of individuals with pathogen loads, i.e. a drift-diffusion equation structured in pathogen load. We are able to prove the existence and uniqueness of the solution to this PDE and to compute explicitly the unique stationary state, which is reached asymptotically. This system is finally generalized by adding an ODE which takes into account the pathogens in the environment seen as a reservoir.

A simple model of the pathobiome dynamics
A simple way to model the pathobiome dynamics in the host gut is to consider two populations: the pathogen and its ecological competitors among the commensal gut bacteria. In order to build the equations, the following mechanisms are considered: the pathogens provoke an inflammatory response of the host, which affects the pathogens as well as the commensal bacteria. However, the inflammation causes an increase in the oxygen level near the gut epithelium, thus providing a competitive advantage to the pathogens. Indeed, the commensal microorganisms are mainly anaerobic, which means that they are highly sensitive to oxygen.
We now detail the model construction step-by-step. We start by assuming an antagonistic relationship between the pathogen and its ecological competitors, i.e. the fraction of microorganisms in the microbiota that share the same ecological niche. Let p and b be the concentrations of pathogens and antagonistic microbes in the host gut, and K > 0 be the carrying capacity of the niche, in the absence of any other limiting factor. The antagonism is modeled by the equation p + b = K (1) Note that p = −b at any time, so that when p increases, b decreases and vice versa, which is consistent with an antagonistic relationship. Note also that (1) allows the elimination of the unknown b in the equations describing the evolution of the system, however for the sake of clarity we will keep both variables until the end of the model description.
The pathogens and the antagonistic bacteria present in the intestine compete to conquer the niche, and the resident bacteria exert a barrier effect against the invasion, meaning that in the absence of any other phenomenon they eliminate the pathogens. Consequently, equation (1) is expanded into the following system of two differential equations: where µ 0 is a positive coefficient accounting for the efficiency of b in the competition.
We now introduce, as mentioned above, a deleterious signal, such as oxygen, resulting from the inflammatory response to the pathogens. Let d be this signal, its evolution in time is modelled as In this equation, the presence of the pathogens directly promotes the host inflammatory response and induces the production of d. This promotion is limited when the number of pathogens become very high; this is expressed through the sigmoid function C 2 p n p n + p n "centered" around the threshold value p ∈ [0, K], with n ∈ N * an exponent controlling the stiffness of the sigmoid. The parameter τ models a time delay between the evolution of the pathogens and the evolution of the signal d.
Two additional phenomena happen. First the deleterious signal d is more harmful to the antagonistic microbiota b than to the pathogens, thus favouring the pathogens in their fight to conquer the niche. Additionally, b undergoes a limitation in the efficiency of its competition with p related to the actual environment in the gut, such as the presence of other competitors of b or the basal immune system response, which we model by a self-inhibition term. These effects are taken into account by modifying the expression of the coefficient µ 0 according to where µ, α and C 1 are positive. The term 1 − α µ b represents the limitation related to the environment, and − C 1 µ d accounts for the influence of the deleterious signal d.
The evolution of b and p is therefore modeled by the following equations in which the antagonistic relationship between b and p is now more complex than a simple competition where b wins and p loses, since asymptotically the introduction of (2) may modify the trade-off between b and p and allow for coexistence. Using (1) to eliminate b and setting A = αK − µ, the overall model is expressed as a system of two equations on p et d: Finally, we assume that the dynamics of d is very fast, that is 0 < τ 1. Therefore, it is natural to assume that d is at quasi-steady state, and to simplify the model into a single equation where C = C 1 C 2 . We then get a monodimensional ODE on the pathogen only, that keeps track of the main drivers of Salmonella infection through its coefficients: 1) the parameter A sums up microbial ecology characteristics of the commensal by balancing its maximal efficiency coefficient µ with the environmental harshness αK, which is also taken into account by the parameter α; 2) the parameter C reflects the virulence of the host inflammatory response (the parameter C 2 describing the amount of deleterious signal produced during inflammation and C 1 its impact on commensals); 3) the host sensitivity to pathogens is shaped by the coefficients p , that determines the amount of pathogen tolerated by the host, and the exponent n that represents the virulence of the inflammation (a higher n induces a sharper inflammatory response to the pathogen presence, beneficiary to the pathogen).
Depending on the values of A, α, C, n and p , this ordinary differential equation can possess from two to five steady states in the interval [0, K], including 0 and K, alternatively stable or unstable. Discarding the trivial situation in which only two steady states (K is stable, 0 is unstable) exist, the stability of the equilibria in the case of 3, 4 or 5 steady states in [0, K] is summarized in Table 1, when A, α and C are positive. As mentioned in the introduction, we are interested in the biologically relevant situation of low and high-shedders, which corresponds to the co-existence of two stable steady states. The equilibria corresponding to 0 (no pathogen in the gut) or K (complete elimination of the antagonistic fraction of the microbiota) are not realistic in a livestock population context, where the contamination by the pathogen is certain and the commensal microbiota resistance is never totally suppressed. So the relevant situation on which we will focus is when there are two stable (p 1 and p 3 ) and three unstable (0, p 2 and K) steady states. The two stable steady states are interpreted as the concentrations of pathogens in low and high-shedders respectively.

A model derived from individual stochastic variability
From the previous ODE model, we obtain now an evolution equation for a whole population of individuals structured by pathogen loads.
Biological systems often exhibit intrinsic stochastic fluctuations in their dynamics. It therefore seems sensible to reformulate the model (3) as a stochastic differential equation (SDE)  Table 1. Stability of the steady states 0 < p 1 < p 2 < p 3 < K when A, α and C are positive.
where P (t) is now a stochastic variable modeling the amount of pathogens in the gut, dB is a gaussian unitary white noise, σ > 0 is the instantaneous standard deviation of the stochastic fluctuations in the model and the function F is defined on R by It can be easily checked that F is uniformly Lipschitz on R, and has linear growth. As σ is constant, we can use standard results from SDE theory (see e.g. [2]) to obtain the so-called forward Kolmogorov equation of (4) the solution u(t, p), defined on L 2 (R + × R), is the probability density function of p(t) at time t, conditionally on the probability density of the initial pathogen population density function u ini . However, this model is not quite satisfactory as brownian stochastic fluctuations make it possible to achieve realizations of p(t) that are negative or greater than K. To avoid this, the SDE (4) has to be transformed into a SDE with reflecting boundary conditions which requires more sophisticated tools from stochastic process theory (see e.g. [3] for a quick introduction). The resulting forward Kolmogorov equation is modified in such a way that the solution space is now L 2 R + , L 2 (0, K) and its formulation is We now consider a very large population such that it can be described by a population density s(t, p), meaning that b a s(t, p)dp is the number of individuals with a pathogen load between a and b at time t. Define s(t, p) = N u(t, p) so that the total size of the population N = K 0 s(t, p)dp is constant. Then the population density s satisfies the following equations: where F is defined at Eq.(5). This is a realistic population model, structured with respect to the pathogen load, in the absence of pathogen transmission.
1.3. Existence of solutions to Eq. (6) and convergence towards a stationary state Let us now prove the existence and uniqueness of the solution to system (6).
Proof. Consider the unbounded operator A with domain such that for all w ∈ D(A), Then −A is a regular Sturm-Liouville operator. It follows from [9] that A generates a strongly continuous semigroup on L 2 (0, K), which proves the proposition.
We now prove that system (6) possesses a unique stationary state.
Proposition 2. The PDE (6) has a unique stationary state defined by where F is defined at Eq. (5).
Proof. The stationary states of model (6) are the solutions of the equation: Thanks to the definition of F and the boundary conditions, this leads to This simple ODE has a unique solution given by formula (7).
Finally, we can show that the solution reaches asymptotically this stationary state. Proof. We first consider solution s(t, ·) of (6) with initial condition in D(A) as defined in proposition 1, and follow the method proposed by Bolley et al. in [7]. We consider the L 2 weighted norm and define the function As the initial condition is in D(A), s(t, ·) is in C 1 ((0, T ), D(A)) and we can differentiate G with respect to t, which leads to Then, integrating by parts, using boundary conditions (6b) and the fact that F (0) = F (K) = 0, see Eq.(5), we find that Then, using the expression of G and G , we obtain that Finally, by a direct application of the differential form of Gronwall's lemma, one can deduce that for all initial condition in D(A), the solution s(t, p) of (6) converges when time goes to infinity towards the steady state s ∞ , in the sense of L 2 ]0, K[ , 1 s∞ , at rate e −σ 2 κt . From the density of D(A) in L 2 (0, K) it follows that for all initial condition in L 2 (0, K), the (mild) solution s(t, p) of (6) also converges when time goes to infinity towards the steady state s ∞ , in the sense of L 2 ]0, K[ , 1 s∞ , at rate e −σ 2 κt .

Generalized model adding transmission through an external reservoir
Our aim in this subsection is to add some exchanges of pathogens between individuals through the environment. The salmonella pathogen can be released in the environment by the infected animals and can contaminate food sources, for example watering troughs. This mechanism creates an external reservoir of pathogens which varies according to the excretion and absorption of salmonella by individuals. Thus, it constitutes a major environmental factor in the spread of pathogens. In order to account for this mechanism, the model (6) is modified into where F is defined at Eq.(5).
Here the new variable r denotes the reservoir of pathogens in the environment. In PDE (8a), the term β in (p)r represents the uptake of pathogens from the environment, while β ex (p) stands for the pathogen excretion by individuals. Similarly, in the ODE (8b), the quantity r(t)s(t, p)β in (p)dp indicates the uptake of pathogens from the environment by individuals with pathogen load p and the last term s(t, p)β ex (p)dp describes the increase of the pathogen reservoir induced by individuals with pathogen load p. We set β in (p) = β in (K − p) and β ex (p) = β ex p, with β in , β ex > 0, so that the contamination from the reservoir decreases with the pathogen load whereas the excretion increases with the pathogen load. The positive parameter γ is the natural decay rate of the pathogen reservoir. Finally, boundary conditions are modified into Robin boundary conditions in order to ensure the conservation of the total population.

Numerical results for models (6) without transmission and model (8) with transmission
Our goal in this section is to display and comment some numerical simulations of systems (6) and (8); we test different initial data and compare the solutions of the two models, with or without transmission.
The numerical scheme is based on upwind scheme for the transport term and the time-implicit centered three point scheme for the diffusion part. The simulations have been performed with Python 3 (using Numpy & Scipy) on a 500 cell mesh grid and with the parameters summarized in table (2).
In this theoretical analysis, in the absence of experimental calibration, the parameters A, α and C are set to values for which the deterministic model (3) has two well separated non trivial stable equilibria corresponding to the last column in Table 1, in order to obtain a bimodal equilibrium distribution in equation (6), describing the coexistence of high-shedder and low-shedder phenotypes in the population. The parameter K has been chosen arbitrarily, and could be rescaled to reflect another carrying capacity. The logistic threshold p has been tuned to trigger an inflammation into the host gut when the pathogen load is close to half the carrying capacity. The diffusion parameter σ was selected in the same order of magnitude than the transport process, meaning that the stochastic effects during the pathogen infection are not negligible and can impact the pathogen dynamics. The parameters β in and β ex provide a small inter-host infection rate.  Table 2. Values of the parameters used to perform the simulations presented in section 2. Figure 1 represents some solutions of the population model (6) structured by pathogen load and without transmission through a reservoir, for the parameters given at Table 2. Three different initial data, shown in dotted blue line in Subfigures 1A, 1B and 1C are considered:

Structured population without transmission
where p 1 ∼ 0.285 and p 3 ∼ 3.714 are the two stable steady states. As expected, for each initial datum, the solution converges towards the theoretical stationary state computed at Eq. (7). This stationary state possesses two maximal values at p = p 1 and p = p 3 , corresponding respectively to the low-shedder and the high-shedder groups and a minimal value at p = p 2 , an unstable steady state of Eq.(3). However, we can observe some discrepancies in the transient behavior, as noticed on the plots of Subfigures 1D, 1E and 1F representing the evolution with time and pathogen load of the density population s. Especially, depending on the initial datum, the steady state can be reached from above (p 3 neighborhood in Fig Figure 2 represents the density population s, solution to the structured in pathogen load population system (8), with transmission through a pathogen reservoir. The initial datum used for this simulation is s ini (p) = 0.2 and r ini = 0.

Structured population with transmission through a pathogen reservoir
The evolution of the population density s(t, p) is plotted in Subfigure 2A, together with levelsets. Numerical experiments show that s(t, p) converges towards a stationary state. The time evolution of the reservoir variable r and of the total pathogen load within the population K 0 ps(t, p)dp is displayed in Subfigure 2B. Due to the pathogen excretion, the reservoir variable first increases until a maximal value of 0.123 at t = 1.11. Then in a second phase, it slightly decreases and tends to stabilize around 0.099. Subfigure 2C shows the comparison between the solution s(t, p) for t = 5, in dashed green line, and t = 10 close to the stationary state (red line). To highlight the effect of the reservoir, we plot s ∞ , i.e. the stationary state without reservoir (r − control, black line), which differs notably from the stationary state with reservoir (red line). It can be observed that low and high shedder clusters are shifted and that their sizes are changed. Indeed, the low shedder cluster is centered around the steady state p 1 0.61 for the model with the reservoir variable, whereas it is centered around the steady state p 1 0.28 without the reservoir variable. In addition, the low shedder cluster is larger in the sense that it contains more individuals when accounting for the reservoir variable. For the high shedder cluster, the effect is the opposite. Indeed, in the model with the reservoir variable, the cluster is smaller (fewer individuals) and shifted towards a smaller value of pathogen load: p 3 3.43 instead of p 3 3.71 without reservoir variable. The increase of the low-shedder average pathogen load may reflect the higher exposure to environmental pathogen, whereas the decrease of high shedder load may come from excretion. Indeed, the degradation of the pathogen in the environment makes the reservoir variable a sink source for intra-host pathogens and pulls the pathogen loads towards lower values.
For this model, stationary states can be computed thanks to some fixed-point technique. However, the proof of the asymptotic convergence towards a stationary state is much more difficult to handle than the proof of Prop.3, due to the coupling with the ODE.
However, we can draw some first conclusions from numerical experiments: for the range of parameters we have tested for β in and β ex , the solution always converges towards a steady state. Moreover, whatever the initial datum s ini chosen among those of Eq. (9) or Dirac masses and whatever r ini , the population density s converges asymptotically towards the same stationary solution.  (6)). Each column corresponds to the different initial data introduced in Eq. (9). Top row: population density is displayed at time t = 0 (blue dotted line), t = 5 (green dashed line) and t = 10 (red dashed line), compared to the stationary state without reservoir (r − control, plain black line). Bottom row: Evolution of the population density, together with level sets, for the three different initial data.
Consequently, we think that the PDE system (8) possesses only one stable stationary solution and that all solutions converge towards this steady state.

Study of various control strategies
Control strategies consist in external actions in order to limit the spread of pathogens within the population. Two types of control will be considered: on the one hand, the cleaning action limits the pathogen reservoir in the environment, by removing a certain amount of pathogens represented by the reservoir variable, and, on the other hand, drugs can be used to cure the population. A combination of these two treatments is also considered and a comparison between these three possibilities is performed.
Let us first describe how to model the two treatment strategies in the PDE system (8).
The cleaning action removes pathogens from the environment. It can be modeled by adding an extra term of the form −C(t, r) in the right hand side of Equation (8b) for the reservoir variable. As a first attempt to describe this mechanism, it is assumed that C is a linear function of r, meaning that the cleaning treatment removes a constant fraction of pathogen reservoir per time unit. Moreover, cleaning might start after a given time and/or might be periodic. Therefore, we assume that C takes the following form: where ρ ∈ R + is the cleaning rate and I C the time interval(s) during which the cleaning is applied.
On the contrary, the treatment with drugs decreases the pathogen load of the individuals. Therefore this mechanism is modeled by adding an extra term −T (t, p)s(t, p) in the drift term of the population mass balance equation (8a). It is assumed that the treatment removes a constant fraction of the pathogen load of an individual per time unit and that it also depends on time. So, we write the term T accounting for drug treatment as where the rate θ ∈ R + models the treatment efficiency and I T the time interval(s) during which the treatment is administrated.
Taking into account cleaning and drug treatment, the system (8) becomes The simulations presented in the following subsections are performed with the parameters of Table 2. We also take ρ = 5 for the cleaning rate and θ = 0.25 for the drug treatment rate.

Cleaning strategy
Let us begin with the study of the cleaning strategy. Figures 3 and 4 represent two simulations of system (10) with cleaning only, i.e. T (t, p) = 0. For the simulation displayed in Figure 3, the cleaning is activated during the whole simulation, namely I C = R + , whereas in the simulation displayed in Figure 4 the cleaning is only activated on the time interval I C = [10, 20].
The results presented in Figure 3 are comparable to the simulation presented in Subsection 2.2. Indeed, in this case, system (10) reduces to system (8) in which the decay rate for the pathogen reservoir γ has been shifted toγ = γ + ρ. As in Section 2.2, the system seems to converge towards a stationary state. However, the features previously observed are enhanced. First, we can notice that the population density s affected by cleaning contains more low shedder individuals and less high shedder individuals than the population density without cleaning. In addition, the reservoir variable takes lower values in the case of cleaning; for example, we observe that at time t = 10, r(10) 0.038 with cleaning and r(10) 0.118 without cleaning. Similarly, the total pathogen load in the whole population takes lower values, for example, at time t = 10, K 0 ps(10, p)dp 2.0 without cleaning and K 0 ps(10, p)dp 1.77 with cleaning. Finally, the pathogen loads corresponding to low and high shedder individuals are changed (p 1 0.42 and p 3 3.39) compared to the model without cleaning (p 1 0.61 and p 3 3.43). Figure 4 displays the results of a simulation where the cleaning is applied only during a range of time I C = [10, 20], which allows to study the dynamic response of the system. In Subfigures 4A and 4B , we can notice three different periods of time. The first one corresponds to the time range t ∈ [0, 10], when the system evolves without cleaning. During this time interval, the system behaves like system (8) (see Section 2.2). The second period corresponds to t ∈ I C = [10, 20], when the cleaning is applied; according to Figure 4B, the pathogen reservoir rapidly drops off to a range of values similar to those found in the simulation where the cleaning is applied from the beginning, see Figure 3B. On the opposite, the dynamic adaptation of the population is slower. Indeed, as it can be observed in Figure 4B, the total pathogen load K 0 ps(t, p)dp, plotted in blue dashed line, also decreases but takes a longer time to stabilize. In the last period, which corresponds to t > 20, cleaning is removed. Here again, the dynamic adaptation of the pathogen reservoir variable occurs very quickly, whereas the adaptation of the population is slower.
As a conclusion, the cleaning action does not change the global dynamics of the system, but leads to smaller reservoir pathogen levels and lower average pathogen load in the population. Moreover, we can observe that the population dynamics response to cleaning takes more time than the pathogen reservoir dynamics response. Furthermore, when cleaning is stopped, the pathogen distribution rapidly turns back to the distribution observed  when no cleaning is applied, showing that cleaning should be applied at a period smaller than this relaxation time scale to be efficient as a control strategy.

Drug treatment strategy
In this section, we now consider the effect of a drug treatment and we neglect the cleaning, i.e. C(t, r) = 0 in Eq. (10). Two different scenarios of treatment are discussed below. Figure 5 represents the first scenario in which the treatment is applied for all time, namely I T = R + . According to Figure 5A displaying the population dynamics and Figure 5B representing the pathogen reservoir, the system seems to converge towards a stationary state. The effect of the treatment can be estimated by comparing this simulation with the results obtained in Section 2.2, that is to say the same system without treatment. In Figure 5B, it can be observed that the pathogen reservoir takes lower values in the case with treatment than without, for example r(10) = 0.053 with treatment and r(10) = 0.10 without. As expected, the total pathogen load in the population is smaller, meaning that the population is less infected by the pathogen when a treatment is administrated. This observation is enforced by comparing the size of the low and high shedder groups in Figure 5C with the ones in Figure 2C. Indeed, with a drug treatment, a large part of the population concentrates around the low shedder pathogen load value p 1 and the high shedder group, i.e. individuals with a pathogen load around p 3 , is smaller. Moreover, the treatment shifts towards lower values the pathogen loads corresponding to low and high shedder individuals, namely (p 1 , p 3 ) = (0.39, 3.09) with the treatment, whereas (p 1 , p 3 ) = (0.61, 3.43) without.  In the simulation represented in Figure 6, drug treatment is only applied during the time interval I T = [10, 20]. As for the cleaning case, this simulation allows to study the dynamic response of the system to the drug treatment and the evolution of the solution in Subfigures 6A and 6B can also be divided into three periods. The first one corresponds to the time range t ∈ [0, 10] with no cleaning. In this first phase, the system behaves like System (8), see Section 2.2 for the corresponding numerical results. The second time period corresponds to t ∈ I C = [10, 20]. During this period, according to Figure 6B, the total pathogen load of the population, plotted in blue dashed line, and the pathogen reservoir decrease exponentially. Remark that, here, both the pathogen reservoir and the total pathogen load within the population decrease with the same exponential rate, unlike what was observed for the cleaning case. In the last period, which corresponds to t > 20, the treatment is removed and, as expected, the pathogen reservoir and the total pathogen load increase. Here again, the dynamic adaptation of the pathogen reservoir and the average pathogen load follow an exponential dynamics with a comparable rate. After the end of the therapy, the system recovers its original dynamics, showing again that the administration periodicity is important for a durable effect.

Combination of cleaning and drug treatment
We now consider the case when cleaning and drug treatment are combined in order to decrease the pathogen load. A simulation of System (10) with both terms C(t, r) and T (t, r) is displayed at Figure 7. In this simulation, cleaning and treatment are activated during the whole time, namely I C = I T = R + . The values for the cleaning rate and the treatment rate are the same as in the previous subsections and, as expected, both effects accumulate. Indeed, by comparing Figure 7B with Figures 3B and 5B, it can be observed that the pathogen reservoir and the  total pathogen load remain lower than in the cleaning or the drug treatment cases. Here again, it seems that the system converges towards a stationary state, according to Figure 7A.

Comparison of the various control strategies
According to the simulations presented in the previous subsections, the cleaning strategy does not allow to decrease efficiently the pathogen load. On the opposite, the drug treatment seems much more efficient and leads to a substantial reduction of the pathogen load in the population. Indeed, cleaning reduces the environmental pathogen load which decreases the positive drift term β in (p)r in Eq. (10a): the impact of cleaning on the pathogen distribution is then driven by the excretion represented by the negative drift β ex (p) which is no longer balanced by the pathogen absorption. The same behavior occurs during drug treatment: the environmental pathogen load is decreased consecutively to the reduction of the high-shedding, and the negative drift β ex also nearly fully applies. However, it is supplemented in that case by the treatment negative drift T (t, p), inducing a stronger impact on the pathogen distribution. Depending on the parameters, the high shedder can be eliminated. However, these conclusions are highly dependent on the chosen parameter values and the modeling of the control strategies is highly simplified. Finally, the combination of cleaning and treatment seems to be, in any case, the most efficient strategy.

Generalization to a compartment model with transfers
In this section, we generalize the model (8) introduced in Sec. 1.4 by considering compartments. In particular, the compartments can represent different cages (or farms), between which there exist exchanges of animal or/and of the pathogens present in the environment, for example pathogens carried by farming tools.
Let us define d ∈ N * as the number of compartments. From now on, will be a vector gathering the populations of the various compartments, structured by pathogen. Similarly, r(t) = r 1 (t), r 2 (t), . . . , r d (t) ∈ R d is now the vector gathering the reservoir variables in all the compartments. Moreover, let A, respectively B, be the transfer matrix in M d (R) representing the exchanges between compartments for the individuals, respectively for the salmonella pathogen reservoirs. Indeed, for 1 i, j d and i = j, the coefficient A i,j 0 (resp. B i,j 0) corresponds to the transfer rate from compartment j to compartment i. Because the total population of animals (resp. pathogens in the reservoir) has to be conserved by the transfers between compartments, the diagonal coefficient A i,i , which corresponds to the total amount leaving the compartment i, satisfies A i,i = − k =i A k,i (resp. B i,i = − k =i B k,i ). With these notations, the Let N (t) = K 0 s(t, p)dp be the vector of the total population in each compartment. By integrating (11a), one can deduce that N satisfies the ODE N (t) = AN (t). Therefore, using Eq.(11d), the population sizes in all compartments are explicitly given by N (t) = exp (tA) N , where N (0) = N = K 0 s ini (p)dp is the vector of the initial total populations in the cages. We also remark that the overall total population remains constant since

Numerical simulations for multiple compartments (cages/farms)
In the following, we consider the case of 4 compartments and we perform some numerical simulations. We will first study the case of transfers of individuals without any transfer of reservoir pathogens, then the case of transfers of reservoir pathogens without any transfer of individuals.

Multiple compartments and transfers of individuals
Here, we investigate the case with no transfers for the reservoir pathogens r, that is to say B = 0, and only transfers of individuals between the compartments. This situation may represent a farm where individuals of the same species will be exchanged between different cages or several farms with exchanges of individuals of the same species between them. Therefore, we make the assumption that the populations have the same characteristics in each compartment, i.e. F , β ex and σ are identical. However, each compartment, cage or farm, may have its own properties and we use therefore different values for the intakes of pathogens from the reservoir β in and for the reservoir decay rate γ according to the compartment. In addition, different initial distributions of the population in terms of pathogen load s ini (p) and various initial conditions of reservoir pathogens r ini are considered. The corresponding values are gathered in Table 3 Note that all the rows of A sum to 0 meaning that there is no loss or gain of individuals induced by transfers. Moreover, we choose s ini such that N = 1 1 1 1 T ∈ Ker(A). Therefore, the repartition of livestock in each compartment remains constant over time. This particular choice has been made to be able to compare the results with the simulations without transfers.  Table 3. Values of the intakes of pathogens from the reservoir β in , reservoir decay rate γ, initial distributions of the population in terms of pathogen load s ini (p) and initial conditions of reservoir pathogens r ini for the 4 different compartments, used in the simulation represented in Figure 8.
The results of the simulation with these parameters are displayed in Figure 8. As it can be observed in Figures 8A and 8B, the system seems to converge towards a stationary state. In Figure 8B the dashed curves represent the pathogen reservoirs in red and the pathogen population loads K 0 ps(t, p)dp in blue as a function of time when the transfers are neglected, whereas the plain curves correspond to the same quantities when the transfers are given by matrix A defined at Eq. (12). Similarly, in Fig. 8C, the black curves correspond to the population distribution s(t = 10, p) as a function of the pathogen load p at time t = 10 when the transfers are neglected; the solutions at different times when transfers are taken into account are displayed, namely s(t = 0, p) in blue, s(t = 5, p) in green, and s(t = 10, p) in red. According to these simulations, the transfers of individuals modify the solution. For example, the pathogen reservoir and the total pathogen load for the population in the first compartment are lower than when transfers are accounted. On the opposite, in the second compartment, the pathogen reservoir and the total pathogen load for the population are higher when the transfers are activated. For the third compartment, the dynamics is almost unchanged. Finally, for the fourth compartment, the values at final time t = 10 for the pathogen reservoir and the total pathogen load of the population are similar but the transient dynamics is slowed down.   (11)). In the different compartment, no particular treatment is applied. Transfer of individuals between compartments is given by the matrix A in Eq. (12), but no transfer between compartment reservoir variables is added (B = 0). Initial data are different in each compartment (see parameter in Table 3 Table 3. We note that consequently, r + control is different in each compartment and that this simulation represents a situation with no inter-compartment transfer (A = 0, B = 0).

Multiple compartments and transfers between pathogen reservoirs
Now, we study the case with no transfers of individuals, that is to say A = 0, and only transfers of reservoir pathogens between the compartments. Transfers of pathogen reservoirs between compartments correspond for example to transfers of pathogens induced by farming activities, for instance by using the same clothes or boots in different buildings on a farm. It can also represent the spread of the pathogen reservoir by environmental factors, e.g. wind, between farms or cages. In this case, all the parameters can be different from one compartment to another. However, in order to compare the effects of the transfers of the pathogen reservoirs with the transfers of individuals, the simulation presented in Figure 9 has been performed with the same parameters as in previous section, see Table 3 According to Figures 9A and 9B, the solution seems to converge in time towards a stationary state. The main effect of transfers of the pathogen reservoirs between compartments can be observed in Figure 9B, which compares the pathogen reservoir size and the total pathogen load in the population over time, in plain curves, with the same simulation without transfers, in dashed curves. With transfer matrix B given at Equation (13), we can notice the effects of the pathogen transfers, for example in compartment 4 or for the time range [0, 2] in compartment 1. However, within the range of parameters used in this analysis, the effect of pathogen reservoir transfers is very weak on the population. Indeed, we can remark that the total pathogen load -see Figure 9B and the stationary states of the pathogen load distribution -see Figure 9C -with and without transfers for the reservoir pathogens are almost identical.

Multiple compartments and spreading of pathogen load
Finally, we study a case with transfers of individuals and reservoir pathogens and we aim at exploring the spreading and contamination of pathogens from the first compartment to the three other compartments. Initial data are therefore given by: and we use the same parameters for all the compartments, that can be found in Table 2. The transfer matrices A and B are given by: such that the first compartment contaminates the second, the second contaminates the third and the third contaminates the fourth. The results are displayed in Fig. 10. As expected, the population densities converge to the same stationary state in the four compartments. With this parameter set, the propagation phenomenon is fast in comparison with the convergence time to stationary state in each compartment. We can observe, in particular in subfigure 10B, that the pathogens reach the next compartment with a time delay, which propagates from one compartment to the other.

Towards a continuous model in space
By considering a large number of compartments and appropriate transfers between them, our model can be seen as a rough approximation of a continuous model in space. For example, the situation where the exchanges between compartments are induced by a diffusion process can be approximated by considering laplacian matrices    (11)). In the different compartments, no particular treatment is applied. Transfer of individuals between compartments is discarded (A = 0) but transfer between compartment reservoir variables is added (B given by Eq. (13)). Initial data are different in each compartment (see parameter in Table 3 Table 3. We note that consequently, r + control is different in each compartment and that this simulation represents a situation with no inter-compartment transfer (A = 0, B = 0).
with periodic boundary conditions for A and B, but with two different diffusion coefficients, namely The parameter values and initial data remain the same as in the two previous subsections and are summarized in Table 3. Here again, we can see at Figure 11 that the solution seems to converge towards a stationary state. Due to the use of diffusion-like transfers, we can observe in Figure 11 that the repartitions of individuals in terms of pathogen loads tend to homogenize. The same kind of homogeneisation can be observed for the reservoir pathogen, its asymptotic values being close in each compartment.  K 0 ps(t, p)dp r(t) with transfers r(t) without transfers K 0 ps(t, p)dp with transfers K 0 ps(t, p)dp without transfers (B) Evolution of the reservoir variable r and total pathogen load within the population in the 4 compartments.   (11)). In the different compartments, no particular treatment is applied. Transfer of individuals and reservoir variables between compartments are considered (A and B given by Eq.(15)). Initial data are given at Eq. (14) and the pathogen load is initially concentrated in the first compartment; the other parameters are independent of the compartments and can be found in Table 2  K 0 ps(t, p)dp r(t) with transfers r(t) without transfers K 0 ps(t, p)dp with transfers K 0 ps(t, p)dp without transfers (B) Evolution of the reservoir variable r and total pathogen load within the population in the 4 compartments.  Figure 11. 4 compartments, diffusive-like transfers: Evolution of the population density in the 4 compartments when contact with an environmental pathogen reservoir is considered (model (11)). In the different compartments, no particular treatment is applied. Diffusive-like transfer of individuals and pathogen reservoir variables is added (A and B given by Eq. (16)). Initial data are different in each compartment (see parameter in Table 3 Table 3. We note that consequently, r + control is different in each compartment and that this simulation represents a situation with no inter-compartment transfer (A = 0, B = 0).