1. Introduction

Consider the mixing of two populations of hosts epidemiologically different with respect to the infection and transmission of a pathogen. What would be the outbreak outcome (e.g., in terms of attack rate) for each host population as a result of mixing in comparison to the situation with zero mixing? To address this question one would need to define what is meant by epidemiologically different and how mixing takes place.

To proceed, let's consider situations where mixing of epidemiologically different populations of hosts occurs. Such situations involve generalist (as opposed to specialist) pathogens capable of infecting multiple hosts and of being transmitted by multiple hosts [33]. Many of such pathogens cause zoonoses such as influenza, sleeping sickness, rabies, Lyme or West Nile, to cite a few [33]. In this paper, we focus on a specific example of a multi-host pathogen, the highly pathogenic avian influenza virus (HPAI) H5N1 -a virus considered as a potential pandemic threat by the scientific community.

The avian influenza virus can infect many hosts: wildfowl and domestic bird species, with occasional spill-over to mammals (including humans); the severity degree of the disease being species dependent: highly lethal (swans, chicken), few deaths (Common Pochards, humans), and asymptomatic (Mallards). Following the re-emergence of the highly pathogenic strain of H5N1 in China 2005 [6,7,28], a series of outbreaks spread throughout Western Europe, including France in 2006 [13,16,20]. The ensuing epizootics showed a need for adapted surveillance programs and a better understanding of the epidemiology of HPAI H5N1 [18]. In this context, this study is part of the French national project for assessing the risk of exposure of domestic birds and poultry farms to avian influenza viruses following introduction by wild birds; although human activities and commercial exchanges are also main sources for introduction of avian influenza [15,17,27,30].

The motivation for this study stems from the 2006 HPAI H5N1 outbreak that took place in France, in the Dombes wetlands. The area is one of the two main routes used by birds migrating across France, and an important stopover, breeding and wintering site for many wild waterfowl species. The outbreak was of minor size and affected mainly wild Anatidae bird species [13,16,20]: Common Pochards (Aythya ferina) and Mute Swans (Cygnus olor). Although the environmental conditions were conducive to the spread of the virus in the Dombes' ecosystem [31,34], it was suggested that the heterogeneity in the response to H5N1 viral infection of different bird species was a possible explanation for the reduced size of the outbreak [13]. Some studies have shown that averaging together different groups of a population, can only lead to a decrease (or no change) observed in the global reproduction number, compared to when no group structure of the population is considered [1]. Ref. [2] pointed out that the variance in the mixing rate between populations can have a substantial effect one the outbreak outcome. Other studies show that for multi-host pathogens, increasing host or species diversity may lead to either reduction or enhancement of the disease risk [12,24]. Therefore, addressing the question posed in the beginning of this section would provide insights and allow advances in the understanding of how avian influenza may spread in such ecosystems.

Our aim in this paper is to use a SIR compartmental model to investigate the effect of host heterogeneity on the disease outbreak in a multi-host population system. More precisely, we study how the outbreak outcome for each constituent population of hosts is affected in a multi-host population system with mixing in comparison with the single-host situation where individual populations are not mixed. The remainder of the paper is as follows. First, the key parameters and response functions characterizing the outbreak outcome are defined and determined for a single-host system in Section 2, and next the defined parameters are used to define the epidemiological heterogeneity in Section 3. Second, Section 4 is devoted to studying how the outbreak outcome in a multi-host population system is changed, due to mixing of epidemiologically heterogeneous hosts, compared to the outbreak outcomes in a single-host situation. Finally, the paper ends with the application of the results in the context of the Dombes area and concluding remarks in Section 5.

---

2. Single-host system

In this section we define the key characteristic parameters of the interacting population-pathogen system and the response function characterizing the outbreak outcome for such a system. To this end, consider a single species or single-host system in which the dynamics of an infection induced by a pathogen can be described within the framework of the compartmental susceptible-infected-recovered (SIR) model ([25]) in which susceptible individuals, $S$, become infected upon contact with infected ones at a rate $\lambda$, infected individuals, $I$, remain infected and are infectious for a mean duration of $1/\alpha$ and they may recover from the infection to become recovered individuals, $R$, with a probability $x$.

At any time $t$, the size of the population is $N=S+I+R$ and the SIR dynamics is described by the system of differential equations given by,

$ \left\{ $$\begin{array}{lcl} \dfrac{{dS}}{{dt}} = -\lambda S \\ [2ex] \dfrac{{dI}}{{dt}} = \lambda S - \alpha I \\ [2ex] \dfrac{{dR}}{{dt}} = x\alpha I \\ [2ex] \end{array}$$ \right. $

(1)

where $\lambda=p\beta I=\beta I/N$ is the force of infection with $p=1/N$ denoting the contact or encounter probability between two individuals, $\beta$ is the rate of infection transmission from an infected to a susceptible and $x$ the fraction of infected individuals that recover. As we will be dealing with situations of short-lived outbreaks, no population renewal by ways of either reproduction or immigration will be considered in this analysis.

In writing Eq.(1) we have used the homogeneously mixing hypothesis and considered that the transmission of infection is frequency-dependent (i.e. the force infection is proportional to the inverse of the population size) like for the true mass-action kinetics [8]. For $x <1$, a proportion $(1-x)$ of the infected population dies, while for $x=1$, the Eq.(1) reduces to the classical SIR model which has been thoroughly studied in the literature [2,8,25].

---

2.1. Characteristic parameters

The above SIR model is characterized by two (non independent) quantities: the generation time $g=1/\alpha$ and the basic reproduction number (i.e., the mean number of secondary cases generated by an index case in a population of susceptible individuals during the entire period of infectiousness) given by [3]¹,

¹The derivation in Ref. [3] goes as follow. Consider a single infected individual applying a constant force of infection $\lambda$ to a naive population of size $N_0$. The probability that any susceptible individual of the population becomes infected at time $t$ is given by, $u(t)=1-{\rm e}^{-\lambda t}$, such that $u(t=0)=0$ at the beginning of the exposure of the naive population to the infectious individual and $u(t\rightarrow\infty)=1$ in the case both the infectious duration and the contact between the infectious individual and the naive population last for a very long time. By the definition of $R_0$, the total expected number of new infected individuals originated from contact with a single infected one is given by the product of the population size times the probability $u(t)$ averaged over the distribution $\psi(t)$ of infectious durations as, $R_0=N_0\int_{0}^{\infty}u(t)\, \psi(t)\, dt$. Now, using an exponential distribution, $\psi(t)=\alpha{\rm e}^{-\alpha t}$, we found, $R_0=\lambda N_0/(\lambda+\alpha)$. The classical expression, $R_0=\lambda N_0/\alpha$, corresponds to the $\lambda\ll\alpha$ limit, and, as it should be, $R_0$ is always smaller than $N_0$ and reaches its maximum $R_0=N_0$ at the $\lambda\rightarrow\infty$ limit. To obtain Eq.(2) we use $\lambda=\beta/N_0$ (as given below Eq.(1) for I=1). Note that Eq.(2) reduces to $R_0=\beta/\alpha$ in the $\alpha N_0\gg\beta$ limit.

$ R_0 = \dfrac{{\beta N_0}}{{\beta + \alpha N_0}}\:, \label{eq:R0} $

(2)

where $N_0=N(t=0)$ is the initial population size. Both $g$ and $R_0$ parameterize the SIR dynamics in Eq.(1) by providing the time and magnitude scales of the infection outbreak. The infection will takeover for $R_0>1$ (major outbreaks) while it will cool down to zero for $R_0\leq 1$ (minor outbreaks). In general, the $R_0$ carries information on the magnitude of the contact-transmission over a period of $g$.

---

2.2. Response function

To define a response function characterizing the outbreak outcome of the SIR model, we consider the following two indicators:

• the reduced persistence or extinction time, $\tau$: starting at $t=0$ the system in an initial state with infected individuals $I>0$, we denote by $t_p$ the time elapsed until the system reaches the state $I=0$ for the first time (representing the persistence of infection in the system until to extinction). We define by the reduced persistence or extinction time the ratio of $t_p$ by the generation time as, $\tau=t_p/g$. The $\tau$ is a random variable with the mean denoted as $\langle\tau\rangle=T$

• the attack rate, $a$: defined as the ratio of the outbreak size (cumulated number of infected individuals) to the initial size $S(0)$ of the susceptible population. The attack rate, $a$, (the fraction of the total number of susceptible who catch the infection during the course of the outbreak) is a random variable given by, $a=1-S(\tau)/S(0)$, where $S(\tau)$ is the susceptible population never infected. The mean value of $a$ will be denoted $\langle a\rangle=A$

To investigate $\tau$ and $a$, we have run SIR stochastic simulations in a population of size, $N_0=2500$ (see Appendix A). Figure 1 illustrates the cumulative distribution functions (with corresponding distributions for $R_0=2$ shown in Fig. 2) of $a$ and $\tau$ for a system with $x = 0$ (all infected individuals die) and $R_0 = 0.5$, $1$ and $2$. Inspection of these figures shows that the distributions of $a$ and $\tau$ are both bimodal with the number and tip values of the modes depending on $R_0$. For the attack rate, the two modes are rather peaked and located around $a=I(0)/N_0$ and $a=1$ with the probability of finding $a$ close to the lower mode (minor epidemics) decreasing with $R_0$. For the reduced extinction time, the values of $\tau$ are smaller for both $R_0 < 1$ and $R_0 > 1$ than for $R_0\sim 1$.

Bearing the distributions of $a$ and $\tau$ in mind, we now deal with the mean attack rate, $A$, and the mean reduced extinction time, $T$, as a function of $R_0$ as displayed in Fig. 3 for $x=0$ (all infected individuals die) and for $x=0.9$ ($90\%$ of infected individuals recover). Clearly, there is a one-to-one functional relationship between $A$ and $T$ indicating that both carry somehow similar information and, therefore, $T$ can be obtained from $A$ and vice versa. Finally, in what follows we will only focus on $A$ for which an approximate expression is derived in Appendix B. Figure. 3 shows that $A$ is null or very small for $R_0 <1$, starts to increase sharply from the threshold $R_0 = 1$ and monotonically increases toward an asymptote of value one for higher $R_0$.

When $x=1$, the contact probability $p$ between two individuals is constant and independent of time, while when $x <1$ the encounter probability between two individuals increases with time as the total population decreases because of deaths. As a result, the mean attack rate increases when $x$ decreases for fixed $R_0$ [e.g., $A(x = 0, R_0)> A(x=0.9, R_0)$] as the force of infection dynamically increases because of the increase in the contact probability when $x$ gets lower.

On the other hand, consider the probability $\omega(I_{th})$ of the occurrence of an outbreak with the total number of infected individuals greater than or equal to a threshold $I_{th}$. Figure 4 shows the results of stochastic simulations of $\omega(I_{th})$ as a function of $R_0$ for different values of $I_{th}$. Interestingly, we find that $\omega(I_{th})$ exhibits a behavior very much alike to the mean attack rate $A$ as a function of $R_0$ and, in addition, $\omega(I_{th})\approx A$ for $I_{th}\geq 10$.

Thus, it follows from what precedes, that the mean attack rate $A$ (and similarly, the probability $\omega(I_{th}\geq 10)$ of an outbreak occurrence) can be considered as the response function characterizing the outbreak outcome for a single host system parameterized by $g$ and $R_0$. Therefore, we define a characteristic response function $F$ as,

$ A =F(R_0, g, x)\:\:;\:\:R_0=F^{-1}[A(g, x)]\:, $

(3)

where $F(\cdots)$ is given in Fig. 3 (or by the approximation in Eq. (16) in Appendix B) and $F^{-1}(\cdots)$ by inverting $F$ from Fig. 3 (or using the approximate Eq. (19) in Appendix B).

---

3. Definition of epidemiological heterogeneity

Within the epidemiological framework as described in the Section 2, a host population interacting with a pathogen can be canonically characterized by two key parameters (or two dimensions): the basic reproduction number, $R_0$, and the generation time, $g$, both describing the tempo and the order of magnitude of an outbreak. Accordingly, two indices $H_R$ and $H_g$ for $R_0$ and $g$, respectively, can be used to characterize the heterogeneity along the two scales of the system such that the overall heterogeneity of the system can be written as, $H^2 = H _R^2 + H _g^2$. For a population of $n$ hosts, with $f_i$ the proportion of the population of host "$i$" ($i=1, 2, \cdots, n$, and $\sum\limits_{i = 1}^n {f_i}=1$), the heterogeneity $H_h$ along a dimension $h$ can be defined as the ratio of the standard deviation to the mean squared:

$ \label{eq:sigma1} H _h = \dfrac{{\sum\limits_{i = 1}^n {f_i h_i^2 } }}{{\left( {\sum\limits_{i = 1}^n {f_i h_i } } \right)^2 }} - 1\:\:;\:\:h_i=R_0, g\:. $

(4)

It follows that a population of $n$ hosts will be considered as epidemiologically heterogenous when the overall $H>0$; the larger $H$ is, the more the system is heterogeneous.

For a single-host population, $H_h=H=0$, and for a two-host system, $n=2$, Eq.(4) can be written as,

$ H _h =y\, \left(\dfrac{z - 1}{zy + 1}\right)^2 {\rm with}{\kern 1pt} y = \dfrac{{f_2}}{f_1} {\rm and}{\kern 1pt} z = \frac{{h_2 }}{{h_1 }} $

(5)

where $h_i$ is either $R_0$ or $g$ for each host. Because of the symmetric relations, $H _h(y, z)=H _h(1/y, 1/z)$ and $H _h(1/y, z)=H _h(y, 1/z)$, the heterogeneity can be calculated anywhere from the $H _h$ in the range $0\leq z\leq1$. $H_h=0$ either in the single-host limit ($y=0$ or $y\rightarrow\infty$) or for $z=1$. Figure 5 shows the reduced one-dimensional heterogeneity, $H_h/y$, as a function of $z$ for different values of $y$. It appears that $H_h/y$ changes a lot as a function of $z$ for fixed $y$ whereas it changes very little as a function of $y$ for fixed $z$.

Note that different demographic fractions $f_i$'s correspond to different $y$, while different $h_i$'s may correspond to the same $z$ [e.g. $z=0.5$ for $(h_1, h_2)=(0.5, 0.25)$, $(1, 0.5)$, $(3, 1.5)$, $\cdots$]. This indicates that, for $y$ fixed, several epidemiologically different situations with identical $z$ may correspond to identical heterogeneity, and $H_h$ does not show any distinction between them.

---

4. $n$-host system

Now, we consider a heterogeneous system (in the sense of Section 3) constituted of $n$ epidemiologically different single-host subsystems interacting with each other by mixing. The question we would like to address here is how the outbreak outcome for each single-host subsystem (characterized as described in the Section 2) will be affected by mixing interactions.

To proceed, consider $n$ single-host subsystems, each of population size $N_i=f_iN$ ($i=1, 2, \cdots, n$), where $N$ is the total population of the system and $f_i$ ($0\leq f_i\leq1$ and such that $\sum\limits_{i=1}^{n}f_i=1$) the population fraction of the host $i$. Within the framework of the SIR model where $S_i$, $I_i$ and $R_i$ denote the number of susceptible, infected-infectious and recovered individuals of the host $i$, with $N_i=S_i+I_i+R_i$, the dynamics of the infection in each single-host subsystem is described by the system of differential equations given by,

$ \left\{ $$\begin{array}{l} \dfrac{dS_i}{dt} = -\lambda_i S_i \\[2ex] \dfrac{dI_i}{dt} = \lambda_i S_i - \alpha_i I_i \\[2ex] \dfrac{dR_i}{dt} = x_i \alpha_i I_i \end{array}$$ \right. $

(6)

where $\lambda_i$, $1/\alpha_i$ and $x_i$ have the same meaning given in Section 2. The mixing interaction between single-host subsystems manifests itself in the force of infection as, $\lambda _i = \sum\limits_{j = 1}^n {p_{ij} \beta _{ij} I_j } $, where $p_{ij}$ is the matrix of contact or encounter probability between two individual hosts $i$ and $j$, and $\beta_{ij}$ is the infection transmission rate in the case of direct contact from an infected host $j$ to a susceptible one $i$.

Assuming a hypothesis of homogeneous mixing of individuals for both within populations of hosts of the same kind (intra) and between host populations of different kind (inter), the elements of the matrix of contact probabilities can be written as,

$ \left\{ $$\begin{array}{l} p_{ii} (t) =\dfrac{1}{N_i(t)}\left[1-\sum\limits_{j=1;j\neq i}^{n}\dfrac{\phi_{ij}N_j(t)}{M_i(t)}\right]\\[2ex] p_{ij} (t) =\dfrac{\phi_{ij}}{M_i(t)}\:\:;\:\:M_i(t)=\sum\limits_{j=1}^{n}\left[1-\delta_{\phi_{ij}, 0}\right]N_j(t) \end{array}$$ \right. $

(7)

where $\delta_{k, l}$ is the Kronecker delta function, equal to $1$ if $k=l$ and to $0$ otherwise, $M_i$ represents the total population entering in contact with hosts $i$ and $\phi_{ij}=\phi_{ji}$ ($0\leq \phi_{ij}\leq 1$ and $\phi_{ii}=1$) is the assortative symmetric matrix that measures the intrinsic affinity to contact between two different kind of hosts $i$ and $j$. Eq. (7) states that intra-host populations are homogeneously mixed and $\phi_{ij}$ measures the degree of homogeneously mixing host populations $i$ and $j$ with $\phi_{ij}=0$ corresponding to zero (inter) contacts between hosts $i$ and $j$ and therefore to non mixed and totally separated populations $i$ and $j$, whereas $\phi_{ij}=1$ corresponding to the case of completely mixed populations where hosts $i$ and $j$ contact each other regardless of their kind. Note that in general $p_{ij}\neq p_{ji}$ and $\displaystyle\sum\limits_{j=1}^{n}p_{ij}N_j=1$.

For the transmission of avian influenza viruses of interest here, we assume that infectious individuals of any kind are efficient sources of virus excretion such that the transmission of the infection to uninfected individuals only depends on the infection susceptibility of the receiver. That is to say that the infection transmission rate $\beta_{i, j} $ from host $j$ to host $i$ only depends on the host $i$, i.e., $\beta_{i, j} = \beta_{i, i} = \beta_i$. In this case, the force of infection can written as,

$ \lambda_i(t) = \left[\dfrac{f_iN_0R_{0, i}}{f_iN_0-R_{0, i}}\right]\, \alpha_i \sum\limits_{j = 1} {p_{ij}(t) I_j(t) }\:\:\mbox{with}\:\:R_{0, i} = \dfrac{{\beta _i f_i N_0}}{{\beta _i + \alpha _i f_i N_0}}\:, $

(8)

where $N_0=N(t=0)$ is the initial population size of the system and $R_{0, i}$ is the intrinsic basic reproduction number of the single-host subsystem $i$ as already defined in Eq. (2).

To go further and for the sake of simplicity, we specialize to the case of $n=2$ subsystems and reformulate the question we are addressing: what would be the outbreak outcome (in terms of attack rate) for each individual subsystem when two epidemiologically heterogeneous ($H \neq 0$) subsystems (each of which is characterized by $R_{0, i}$ and $g_i$) are epidemiologically interacting ($\phi_{ij} \neq 0$)?

For the mixing between $n=2$ epidemiologically different single-host subsystems, the population fractions are such that $f_1+f_2=1$ and the assortative matrix reduces to $\phi_{ij}=\phi$. The two-host population SIR system thus involves seven independent parameters: $N$, $f_1$, $R_{0, 1}$, $R_{0, 2}$, $\alpha_1$, $\alpha_2$ and $\phi$. The outbreak outcome can be analyzed at two levels: the whole system level and the subsystems level of each constituent.

---

4.1. Outbreak outcome at the whole system level: Global reproduction number, ${\mathcal R}_0$

General considerations on the outbreak outcome can be drawn from the ${\mathcal R}_0$ of the entire system. The disease will invade the multi-host population when ${\mathcal R}_0>1$ but it will stutter to die for ${\mathcal R}_0 <1$ [9,10,21]. The threshold ${\mathcal R}_0$ for a multi-host system can be determined using the next generation matrix (NGM) approach [11]. In the case of a two-host system, the NGM, $\mathbf{K}$, is a $2\times 2$ matrix with the elements given by,

$ \label{eq:IA3} \left\{ { $$\begin{array}{lcl} {K_{1, 1} = \left( {\dfrac{{R_{0, 1} f_1N_0}}{{f_1N_0 - R_{0, 1} }}} \right)\left[{1-\phi f_2} \right]} ; {K_{1, 2} = \left( {\dfrac{{R_{0, 1} f_1N_0}}{{f_1N_0 - R_{0, 1} }}} \right)\left(\dfrac{{\alpha _1 }}{{\alpha _2 }}\right)\phi f_1} \\\\[2ex] {K_{2, 1} = \left( {\dfrac{{R_{0, 2} f_2N_0}}{{f_2N_0 - R_{0, 2} }}} \right)\left(\dfrac{{\alpha _2 }}{{\alpha _1 }}\right)\phi f_2} ; {K_{2, 2} = \left( {\dfrac{{R_{0, 2} f_2N_0}}{{f_2N_0 - R_{0, 2} }}} \right)(1 - \phi f_1)} \\ \end{array}$$ } \right. $

(9)

In this approach, ${\mathcal R}_0$ corresponds to the dominant eigenvalue of $\mathbf{K}$, given by,

$ \label{eq:eigenvalue} {\mathcal R}_0= \frac{1}{2}\left[{K_{2, 2} + K_{1, 1} + \sqrt {\left( {K_{2, 2}-K_{1, 1} } \right)^2 + 4\left( {K_{2, 1} K_{1, 2} } \right)} } \right]\:. $

(10)

Because of the term $K_{2, 1}\times K_{1, 2}$ in Eq.(10), ${\mathcal R}_0$ is independent of $\alpha_1$ and $\alpha_2$ indicating that heterogeneity $H_g$ alone do not impact the outbreak. A subsequent sensitivity analysis was conducted using the extended Fourier amplitude sensitivity test (FAST) [32] to study the effects of the parameters $N$, $f_1$, $R_{0, 1}$, $R_{0, 2}$ and $\phi$ on ${\mathcal R}_0$. Figure 6 summarizes the results of the sensitivity analysis. It appears that the parameters having both main and total effects (main effect plus interaction with other parameters) on ${\mathcal R}_0$ are (ranked from the least to the most important): $N$, $\phi$, $f_1$, $R_{0, 1}$ and $R_{0, 2}$; indicating that the outbreak outcome is sensitive to host heterogeneity (characterized at least by different $R_{0, 1}$ and $R_{0, 2}$, i.e., $H_R\neq 0$). Inspection of Eq.(10) indicates that:

• For a fixed nonzero heterogeneity $H_R$, the ${\mathcal R}_0$ always decreases as a function of $\phi$ from ${\mathcal R}_0 = max(R_{0, 1}, R_{0, 2})$ (by definition of the dominant eigenvalue) at zero mixing $\phi=0$ to ${\mathcal R}_0 = R_m$ at complete mixing $\phi=1$, where $R_m$ is the population weighted average of intrinsic reproduction numbers given by:

$ R_m= \left( {\dfrac{{ f_1N_0}}{{f_1N_0 - R_{0, 1} }}} \right)R_{0, 1}f_1 + \left( {\dfrac{{f_2N_0}}{{f_2N_0 - R_{0, 2} }}} \right)R_{0, 2}f_2 \:. $

(11)

The decreasing of ${\mathcal R}_0$ with $\phi$ is due to the crowding or herd immunity effect. Note that at exactly zero mixing the infection remains confined in the host population initially infected, i.e., ${\mathcal R}_0 = R_{0, i}\neq \mbox{max}(R_{0, 1}, R_{0, 2})$ where $i$ corresponds to the initially infected host population.

• For a fixed nonzero mixing $\phi$, the effect of $H_R$ on ${\mathcal R}_0$ is not straightforward because ${\mathcal R}_0$ is neither an explicit nor an implicit function of $H_R$. Since $H_R$ is a function of demography $y$ and reproduction $z$ (see Section 3), two scenarios can be drawn:

- for any fixed ratio of reproductive numbers $z$, the ${\mathcal R}_0$ decreases as a function of $y$ (i.e., $H_R$) due to the crowding effect

- for fixed demography $y$, the ${\mathcal R}_0$ increases with $max(R_{0, 1}, R_{0, 2})$ (rather than the ratio $z$).

The ${\mathcal R}_0$-level curves in Fig. 7 illustrate the effects of $H_R$ and $\phi$ on the variety of behaviors of ${\mathcal R}_0$ as a function of $R_{0, 1}$ and $R_{0, 2}$.

---

4.2. Outbreak outcomes at the subsystem level: Equivalent basic reproduction number, $R_{eqv, i}$.

To investigate the effects of mixing on individual outbreak outcomes at the level of each subsystem, we have run SIR stochastic simulations in a two-host system (see Appendix A) with a total population of size, $N_0=5000$ and heterogeneities, $H_R\neq 0$ and $H_g = 0$. For the purpose of the investigation in the context of avian influenza, we have set the parameters $\alpha_1=\alpha_2=0.2$, $x_1 = 0$, $x_2 = 0.9$ while the others $\phi$, $f_1$, $R_{0, 1}$ and $R_{0, 2}$ are allowed to vary. Values of $\alpha$'s and $x$'s are chosen to correspond to mean durations of virus excretion and recovering probabilities for swans (highly susceptible species) and ducks (mildly susceptible species) infected with HPAI H5N1 [4,5,22,23,26,29].

Figure 8 illustrates the cumulative distribution (cdf) of the attack rates for each host in the system and for the whole system. The cdf of the whole system is broad and close to that of the most abundant population host $1$ ($f_1>f_2$). The distribution of attack rate is almost bimodal for host $1$ (similar to a single-host system) while it is much broader for the host $2$ as a result of mixing. In what follows, we focus on the mean attack rates.

Because of mixing, the mean attack rate $A_i$ for each host $i$ involves two contributions: $A_i=A_{ii} + A_{ij, i\neq j}$, where $A_{ij}$ is the mean attack rate in the host population "$i$" caused by infections transmitted by the host population "$j$", with $A_{ii}<A_{0, i}$ where $A_{0, i}$ is the mean attack rate at zero mixing or for a single-host system. To assess to which extent mixing and heterogeneity affect the outbreak response for each host population, we consider the ratio,

$ \eta_i=\dfrac{F_i^{-1}(A_i)}{F_i^{-1}(A_{0, i})}=\dfrac{R_{eqv, i}}{R_{0, i}}\:, $

(12)

where we have used the relation in Eq.(3) (see Section 2) to define the equivalent basic reproduction number as, $R_{eqv, i} = F_i^{-1}(A_i)$, where $F_i ^{-1}$ is obtained from the numerical inversion of $F_i $ given in Fig. 3. No effect (i.e., $A_i=A_{0, i}$) corresponds to $\eta_i=1$, while $\eta_i>1$ and $\eta_i <1$ indicates amplification and dilution effects for host "$i$", respectively. The amplification (dilution) effect occurs when the outbreak response for species "$i$" in the multi-host system is greater (smaller) than the one in the absence of population mixing [24].

Several combinations of $R_{0, 1}$ and $R_{0, 2}$ with varying $f_1$ and $\phi$ were explored using stochastic simulations (see Appendix A). Results are summarized in Table 1, showing the different kinds of outbreak responses observed when mixing two host populations, for which the infection transmission between individuals depends on the infection susceptibility of the receiver:

• three kinds of behaviors for each host population are possible depending on the mixing and heterogeneity parameters: dilution, no effect or amplification behaviors. As shown in Table 1, the interaction between two heterogenous hosts, with at least a $R_{0, i} > 1$, leads to a dilution effect for the hosts with higher $R_{0, i}$ and to an amplification effect for the hosts with lower $R_{0, i}$. In terms of metaphor, this is reminiscent to the thermalization effect when mixing two miscible liquids at different temperatures. However, this effect does not hold when $R_{0, i} < 1$ for the two hosts where dilution effects are observed for both hosts regardless their relative values of $R_{0, i}$. The dilution effect is due to the crowding or herd immunity effect while the amplification stems from the addition of cross infected cases ($A_{ij}$) between host subsystems, and no effect results from either zero impact or compensation of dilution and amplification effects. Therefore, the whole system characterized by ${\mathcal R}_0$ consists of the coexistence of two-phase behaviors, i.e., subpopulation of hosts undergoing a dilution effect while the other one an amplification effect.

• the extent to which a subsystem undergoes dilution or amplification is a function of demographic and mixing parameters with a possible transition from dilution via no effect to the amplification behaviors (and vice versa), when varying the individual $R_{0, i}$.

• as the proportion of recovered $x_i$ did not appear explicitly in $R_{0, i}$, we purposely did not include $x_i$ in the heterogeneity indices. Table 1 shows that the overall heterogeneity does matter in the outbreak outcome indicating hence the need to incorporate $x$ into $H$.

---

4.3. Overall effect of heterogeneity

Figures 9 and 10 illustrate some of the situations presented in Table 1. Figure 9 shows the coexistence of two-phase behaviors (dilution effect for a subpopulation and amplification effect for the other one), where the ${\mathcal R}_0$ of the whole system increases with $\mbox{max}(R_{0, 1}, R_{0, 2})$ and decreases with $\phi$, whereas outbreak responses of host $1$ (dilution or amplification) increase with $\phi$. In Fig 10, we explore the effect of initial conditions ($I_{i}(0)$, i.e. in which host population the infection starts) not shown in Table 1. Clearly ${\mathcal R}_0$ is not sensitive to initial conditions and decreases with $\phi$ as expected, whereas the extent of the two coexisting behaviors (dilution effect for a subpopulation and amplification effect for the other one) very much depends on the initial conditions for low $\phi$. The extent of dilution and amplification effects for the same system with same parameters may be different at low $\phi$ depending on initial conditions and become identical at higher $\phi$.

---

5. Concluding remarks

The aims of this work were to define the epidemiological host heterogeneity and investigate the effect of host heterogeneity on the disease outbreak outcomes for each host in a multi-host population system, given prior knowledge of the disease epidemiology for each host population in the zero mixing situation. In other words, what is the impact of a multi-host system on the outbreak response of individual host populations involved?

We have shown that a single-host system can be canonically parametrized using two quantities, the basic reproductive number $R_0$ and the generation time $g$, and characterized by an epidemic or outbreak response function $F(R_0, g, x)$ (like in Fig. 3) describing how a host population responds (in terms of attack rate, persistence time) to an infection introduction. To deal with a heterogeneous multi-host system involving epidemiologically different single-host populations, two ingredients must be considered:

• Heterogeneity index $H$: involving two dimensions, the generation time and the basic reproduction number, $H$ measures to which extent host populations are different in terms of $R_0$ and $g$ combined with their demographic weights $f$. A homogeneous system correspond to $H=0$ while $H>0$ corresponds to epidemiologically different populations, having different epidemic response functions $F(R_0, g, x)$.

• Interaction matrix: which takes into account both epidemic and demographic characteristics to structure how different hosts interact with each other. By interactions we mean that hosts have an epidemic and a demographic role in the transmission and spreading of the infection. For the two-host case presented in this analysis, the control parameter for the interaction matrix reduces to a single assortative mixing index $\phi$ that measures the degree of homogeneously mixing two kinds of host populations.

As minimal definition and necessary conditions, we state that the epidemiological host heterogeneity occurs in a system of epidemiologically interacting populations where each host population is characterized by a different epidemic response function. There is no host heterogeneity in the absence of interactions between populations or when interacting populations have all identical epidemic response functions.

Regarding the impacts of host heterogeneity on the outbreak outcomes, we found that they are twofold in the case of the infection transmission depending on the receiver infection susceptibility: i) -outbreak dampening, i.e., the outbreak in the heterogeneous multi-host system is always smaller than the summation of outbreaks for individual subsystems taken separately, and ii) -as summarized in Table 1, three kinds of outbreak outcomes are possible for the individual subsystem depending on the mixing and heterogeneity parameters: dilution, no effect or amplification behaviors where the outbreak responses in the multi-host system are lower, similar or higher than in the single host system, respectively, with the magnitude depending both on $H$ and $\phi$.

Previous works, [14], have shown that, in the case of preferential mixing, like in this study (though with a different mixing pattern), the disease can invade the population when any subgroup is self-sufficient for the disease transmission (i.e., $R_{0, i} >1$). In addition, increasing the intra-group mixing rate increases the probability for the disease to invade the population. This is consistent with our finding that the global ${\mathcal R}_0$ decreases when $\phi$ increases (i.e., decreasing the intra-group mixing rate). However, in our model, having an individual $R_{0, i}>1$ is not sufficient to ensure ${\mathcal R}_0 >1$. Indeed, for $\phi=1$, the global ${\mathcal R}_0$ is given by Eq. 11, which can result to ${\mathcal R}_0 <1$ for sets of $f_1$, $f_2$, $R_{0, 1}$ and $R_{0, 2}$. This difference stems from differences in the mixing structure between our model and that described in [14].

The previous works were largely focused on the impacts that heterogeneity may have on the global ${\mathcal R}_0$, or the ability of a disease to invade a population consisting of different subgroups, as opposed to a homogenous population. Even though we address this issue in our paper as well, we highlight the effects that a mixing of heterogenous sub-populations has on the ability of the disease to invade each sub-population, compared to a situation were no mixing is considered.

The situation of the HPAI H5N1 outbreak in mid-February 2006 in the Dombes, France, can be analyzed within the framework of the afore outlined approach. As mentioned in the Introduction section, although the environmental conditions were conducive to the spread of the virus in the Dombes' ecosystem [31,34], the outbreak was of minor size, mainly affecting Common Pochards (Aythya ferina) and Mute Swans (Cygnus olor) [13,16,20]. It was suggested that the host heterogeneity in the response to H5N1 viral infection of different bird species was a possible explanation for the reduced size of the outbreak [13].

During the outbreak period, the situation in the Dombes was that Swans, Common Pochards and Mallards were found well mixed with a census of $600$, $15000$, and $7500$, respectively [20]. Given that Swans are highly susceptible to influenza virus infection with a short mean death time, and a high viral excretion level, Common Pochards are less susceptible than Swans to influenza virus infection, with low mortality rate, and Mallards have a low susceptibility to influenza virus infection with no associated mortality [13], we may infer that $R_{0, swan}>R_{0, pochard}>R_{0, mallard}$ with likely $R_{0, swan}>1$. In addition, analysis of migration patterns indicated that as the swan population migrated to the Dombes area about two months before the outbreak onset, the disease was therefore likely introduced in the area by ducks (Pochards) following their massive arrival in early February as a result of a cold spell [19]. Based on all of this, the 2006 outbreak in the Dombes can be described within the framework of epidemiological host heterogeneity with $R_{0, duck} <1$ (combining Common Pochards and Mallards) and $R_{0, swan}>1$ just like llustrated in Fig 10 (with Swans = 1 and ducks = 2 and initial conditions $I_1(0)=0$ and $I_2(0)=1$). In this figure, Swans (host 1) undergo a dilution effect, while ducks (host 2) an amplification effect. Given the small numbers of both dead Swans and infected and dead ducks during the outbreak, one may suggest that $i$) -the duck $R_{0, duck}$ was substantially smaller than one, $ii$) -the mixing $\phi$ between Swans and ducks was large perhaps close to one, and $iii$) -even speculate that the epidemic threshold transition from $R_{0, swan}>1$ to $R_{eqv, swan} <1$ occurred for swans.

To conclude, we have depicted a framework for defining the epidemiological host heterogeneity and assessing its impacts on outbreak outcomes in terms of epidemic response functions for host populations in interaction. The approach was illustrated for the case of frequency-dependent direct transmission where the infection transmission depends on the receiver infection susceptibility, (i.e., $\beta_{i, j} $ from $j$ to host $i$ only depends on the host $i$, i.e., $\beta_{i, j} = \beta_{i, i} $) and the two-host system was used as the minimal multi-host system. This work can be extended in several other directions: generalization to $n>2$ hosts systems, use of a general transmission matrix $\beta_{i, j} $, and including spatial heterogeneities.

---

Appendix A. Stochastic simulations of SIR

Stochastic simulations for the SIR model were generated using the stochastic discrete time version of the system of equations in Eq.(6), in which $S_i$, $I_i$ and $R_i$ for the each host $i$ are random variables with the transitions:

$ \left\{ $$\begin{array}{lcll} (S_i, I_i, R_i) & \rightarrow & (S_i -1, I_i +1, R_i) & \mbox{at rate } \lambda_i(t) S_i \\ [2ex] (S_i, I_i, R_i) & \rightarrow & (S_i, I_i -1, R_i +1) & \mbox{at rate } \alpha_i I_i \mbox{ with probability } x_i\\ [2ex] (S_i, I_i, R_i) & \rightarrow & (S_i, I_i -1, R_i) & \mbox{at rate } \alpha_i I_i \mbox{ with probability } 1-x_i \end{array}$$ \right. $

(13)

describing the transition from susceptible to infected following a Poisson process of parameter $\lambda_i(t)$, the sojourn times in infected-infectious state by an exponential distribution of mean $1/\alpha_i$, and the probability for infected to recover by $x_i$. To avoid uncontrolled changes in $\lambda_i(t)$, the step "$(S_i, I_i, R_i) \rightarrow (S_i -1, I_i +1, R_i)$" must be performed in parallel for all hosts prior to others processes. Starting from the initial conditions, $(S_i(0), I_i(0), R_i(0))$, at any time $t$ the population for a given stochastic trajectory is given by the random vectors, $(S_i(t), I_i(t), R_i(t))$. All calculations are implemented in Matlab software, release 7.0

• Single-host system: The subscript $i$ can be dropped and the $\lambda(t)$ is expressed in terms of $R_0$ using Eq.(2) as,

$ \lambda(t) = p\beta I=\left[\dfrac{N_0R_0}{N_0-R_0}\right] \alpha\times \dfrac{I(t)}{N(t)}, $

(14)

where $N(t)=S(t)+I(t)+R(t)$. In addition, the time is scaled by $1/\alpha$. A total of $4\times 10^4$ stochastic simulations in a population of size, $N_0=2500$, with the initial conditions, $(S(0), I(0), R(0))=(N_0-I(0), 1, 0)$.

• Two-hosts system: $\lambda_i(t)$ is given by Eq.(8),

$ \lambda_i(t) = \left[\dfrac{f_iN_0R_{0, i}}{f_iN_0-R_{0, i}}\right] \alpha_i \sum\limits_{j = 1} {p_{ij}(t) I_j(t) }, $

(15)

where $p_{ij}(t)$ are the time-dependent contact probabilities. A total of $4\times 10^4$ stochastic simulations in a population of size, $N_0=5000$.

---

Appendix B. Approximate expression of the mean attack rate

When all infected individuals recover from infection, i.e., $x=1$, an equation for $A$ can be derived from Eq.(1) as ([25]),

$ A = 1-\exp\left\{-\left(\dfrac{R_0}{N_0-R_0}\right) \left[I(0)+AS(0)\right]\right\} . $

(16)

For $x <1$ there is no simple way for deriving an expression for $A$. However, an approximation of $A$ can be derived using the two modes, $a=I(0)/N_0$ and $a=1$, of $a$ as follows,

$ A= \left(\frac{I(0)}{N_0}\right)\times u + 1\times (1-u), $

(17)

where $u$ can be regarded as the probability of minor epidemics. We found by numerical analysis that $u$ can be described by,

$ u= \tanh\left(c \times e^{ - bR_0 }\right) $

(18)

where the constants $b$ and $c$ depend on $g$, $x$ and $N_0$. Figure 11 shows the comparisons between simulation results and analytical expressions for $u$ and $A$ given by Eqs. (17) and (18), respectively, with $c = 10.375$ and $b = 2.123$ for $x=0$. Therefore, the expression given in Eq. (17) is considered as an approximate of the characteristic response function $F$, and the inverted function $F^{-1}$ is given by,

$ R_0 =F^{-1}(A)= - \frac{1}{b}\ln \left\{ - \frac{1}{2c}\ln \left[\frac{I(0)-A N_0}{I(0)-(2-A) N_0}\right]\right\} . $

(19)

---

Acknowledgments

AM is a PhD student supported by a grant from the Ministry of Education and Research of France through the Ecole Doctorale Ingénierie pour la Santé, la Cognition et l'Environnement (EDISCE) of Grenoble Alpes University.We are grateful to M. Artois for fruitful discussions. This work has benefited from the support of the Ministry of Agriculture and fisheries under the Project Cas DAR 7074.

---