key: cord-0776367-zbfd03f6
authors: Borovkov, K.; Day, R.; Rice, T.
title: High host density favors greater virulence: a model of parasite–host dynamics based on multi-type branching processes
date: 2012-03-30
journal: J Math Biol
DOI: 10.1007/s00285-012-0526-9
sha: f9bab971b7a0cf9d130f86dcb439981f280e2276
doc_id: 776367
cord_uid: zbfd03f6

We use a multitype continuous time Markov branching process model to describe the dynamics of the spread of parasites of two types that can mutate into each other in a common host population. While most mathematical models for the virulence of infectious diseases focus on the interplay between the dynamics of host populations and the optimal characteristics for the success of the pathogen, our model focuses on how pathogen characteristics may change at the start of an epidemic, before the density of susceptible hosts decline. We envisage animal husbandry situations where hosts are at very high density and epidemics are curtailed before host densities are much reduced. The use of three pathogen characteristics: lethality, transmissibility and mutability allows us to investigate the interplay of these in relation to host density. We provide some numerical illustrations and discuss the effects of the size of the enclosure containing the host population on the encounter rate in our model that plays the key role in determining what pathogen type will eventually prevail. We also present a multistage extension of the model to situations where there are several populations and parasites can be transmitted from one of them to another. We conclude that animal husbandry situations with high stock densities will lead to very rapid increases in virulence, where virulent strains are either more transmissible or favoured by mutation. Further the process is affected by the nature of the farm enclosures.

Mathematics Subject Classification Primary 60J85; Secondary 90D30 · 60J80 · 92B99

There is an extensive literature on mathematical modeling of infectious diseases. Both the rate of spread and the virulence of pathogens are important, and both the dynamics of disease spread and the evolution of virulence in parasite-host systems have attracted much attention. The first mathematical model of disease dynamics apparently goes back to the 1760 paper by Bernoulli (1760) (see also pp. 2-6 in Daley and Gani 1999) . A major step in further development of deterministic modelling was the "threshold theorem" (Kermack and McKendrick 1927) showing that the initial frequency of susceptibles must exceed a certain threshold value to give rise to an epidemic. The first stochastic epidemiological models appeared in the late 1920s (McKendrick 1926) , applying the "law of mass action" suggested in Bernoulli (1760) to probabilistic description of the epidemics. Stochastic models became much more popular after Bartlett (1949) formulated a model for the general stochastic epidemic by analogy with the deterministic model from Kermack and McKendrick (1927) . Anderson and May (1982) introduced the first model of evolutionary change in pathogen virulence in 1982, which explained the decline in virulence of myxomatosis, introduced to control Australian rabbits, in terms of optimal transmission of the parasite between hosts. Subsequent work, notably Ewald (1983) , Massad (1987) , Knolle (1989) , Lenski and May (1994) , Gandon et al. (2001) , has expanded on this classic paper, to cover the effects of a vector of the pathogen, vaccination of hosts, and various other issues. For more extensive surveys of the field and literature reviews, we refer the reader to the monographic literature (Anderson and May 1991; Becker 1989; Daley and Gani 1999; Andersson and Britton 2000; Mode and Sleeman 2000) and also to papers (Allen 2008; Dietz and Schenke 1985; Hethcote 1994; Ebert 1999; Day and Proulx 2004) .

In natural systems, the density of host organisms declines when a pathogen produces an epidemic. Almost all the models referred to above have considered such natural populations of hosts, where the density of susceptible hosts is reduced by the pathogen. They show that as host density is reduced, less virulent strains of pathogen become more viable, such that infected hosts live longer, and thus tend to pass the pathogen to more susceptible new hosts, maintaining the epidemic. While hosts will have longer generation times than their pathogens, they too, will evolve to be more resistant, further reducing the virulence of the pathogen for these hosts. On this basis one might expect pathogens to be highly virulent only when they have recently switched hosts, such as the SARS coronavirus that was probably transferred to humans from animals such as Himalayan civet cats (Guan et al. 2003) or when the plague bacterium, Yersinia pestis, has been transferred to humans from rodents (Yersin 1894). Various reasons, such as low host resistance, have been suggested to explain high initial virulence of transferred pathogens (Bolker et al. 2010) .

Most models have used a game-theoretic approach (see e.g. Anderson and May 1982; Day 2002; Frank 1996) and focused on the equilibrium level of both host density and pathogen virulence in natural host populations infected by a pathogen. But a number of papers have noted that the initial epidemic phase of pathogen infection may select for high virulence (Lenski and May 1994; Frank 1996; Day and Proulx 2004; Bull and Ebert 2008; Bolker et al. 2010) . In particular, Day and Gandon (2006) , Day and Proulx (2004) present a general theory of evolution of virulence that is capable of predicting both the short-and long-term evolution of virulence.

In animal husbandry as opposed to natural populations, epidemics of new pathogens are not left to reduce the densities of the stock. Instead, antibiotics are applied, hosts are destroyed, or other measures are used in an attempt to eliminate the pathogen, once it is recognised that a serious pathogen epidemic is spreading. If such pathogens persist, they must do so through resistant stages and some mechanism of transfer to new enclosures of domestic stock or crops, where the initial epidemic stage is repeated. The initial stages of each epidemic will typically occur at very high host density. This is the situation we wish to model.

The main objective of the present paper is to provide quantitative and also qualitative [cf. the discussion of the above-mentioned threshold theorem and other findings from Kermack and McKendrick (1927) on pp. 11 and 29 in Daley and Gani (1999) ] insights into the dynamics of pathogen populations during the initial stages of epidemics, assuming that in this period the virulence of the pathogens can change due to mutation. More specifically, we are interested in analyzing the effects that host density may have on the virulence of the pathogen.

In animal husbandry, stock are commonly maintained at very high densities. Thus the question arises as to whether high host densities are likely to lead to selection for more virulent strains of pathogens during the initial stages of epidemics. This question has been investigated only recently (Day and Proulx 2004; Bolker et al. 2010 ). Both these studies use models of the epidemiological dynamics of natural host populations and the evolution of the pathogen, incorporating an assumed tradeoff between virulence (the extra mortality due to infection) and transmissibility (the probability of transmission from an infected to a susceptible host in contact with it). We propose a novel approach to this question.

The mathematical models we discuss in the present paper incorporate a continuous time Markovian multitype branching process. Branching processes first appeared in epidemiological context in the mid-1950s: one can mention here the paper by Bharucha-Reid (1958) and also refer to Bartlett (1955) , Kendall (1956) , Whittle (1955) . For further literature review, we again refer the reader to the above-listed monographic and survey literature, while some examples of applications of multitype branching processes in epidemiology can be found in Becker and Marschner (1990) , Meester et al. (2002) , Yates et al. (2006 ), Heinzmann (2009 ), Alexander and Day (2010 .

Branching processes are relatively simple and well-understood mathematical models; for treatments of the theory of branching processes, see e.g. Jagers (1975) , Athreya and Ney (1972) , a recent addition to the monographic literature on branching processes prepared specifically for biological audience being Haccou et al. (2005) . At the same time, the processes were shown to be good approximations to the general stochastic epidemic models at the initial stage (when the total population size is large and the initial number of infectives is small) and also at the final stages of the epidemics (see e.g. Ludwig 1975 , Ball and Donelly 1995 , Hethcote 1994 . As our main objective is not to model the detailed dynamics and describe the whole history of the epidemic itself, but rather to suggest a model for changes in the pathogen's population composition depending on the density of the host population, we will restrict ourselves to working with the simpler branching process models.

One can also envisage an extension of our simple analysis to the general stochastic models as the same mechanism will certainly work for the latter as well. However, such extensions will be much less tractable analytically and may lead to no closed-form answers.

Moreover, the most critical stage of an epidemic is the initial one, when it is basically determined if there will be a large-scale event or the epidemic will die out. And as branching processes are good approximations to the general stochastic epidemic models at the initial stage, the threshold analysis aimed to determine if the "basic reproductive number" (defined roughly as the expected number of secondary cases produced by one infected individual) is greater than one (which implies the danger of an outburst or persistence of endemic levels) can be carried out using those models (see e.g. Chapters 6 and 8 in Mode and Sleeman 2000) .

In our analysis, we will look at supercritical two-type branching processes (so that the basic reproductive number will be greater than one: we are interested in what happens when there exists a threat of epidemics) and then look at the behaviour of the ratio of the sizes of the subpopulations in the process (representing two versions of the pathogen, that can mutate into each other). This quantity can be used to determine which of the two types will become dominant in the population over time.

As we are interested in stock or crop systems of high host density not regulated by the parasite, cf. Lenski and May (1994) , we assume natural mortality of the hosts is zero and that host density is effectively not changed by the infection during the period of interest. We define the lethality of the pathogen as the mean survival time of an infected host [which corresponds to the usual definition of virulence as the extra mortality due to infection (Day 2002; Day and Proulx 2004; Bolker et al. 2010) ], and the transmissibility as the probability of an infected host infecting a susceptible one if a contact occurs [which corresponds to the transmission rate per susceptible host in many published models Bolker et al. (2010) ]. For recent discussions of the interplay between these two characteristics, see Girvan1 et al. (2002) , Lipsitch and Moxon (1997) . The theory of a tradeoff between these postulates that high virulence (which involves higher host exploitation) requires high transmission, but this theory is still debated (Alizon et al. 2009; Bolker et al. 2010) .

In the paper, we use our approach to model the dynamics of a host/parasite population where parasites can be of one of two types that can differ in their lethality and transmissibility. We do not assume there is any necessary relation between these, as in most models (see Bolker et al. 2010) . The underlying simple continuous time Markov model of a two-type branching process is presented in Sect. 2. The analysis of the model and derivation of the dynamics for the mean functions are given in Sect. 3 and used in Sect. 4 to establish the eventual composition of the pathogen population. It turns out that our model does show changes in the overall parasite's lethality in response to increased density of the hosts. This is not dependent on any relationship between the lethality and transmissibility of a pathogen strain, as assumed in many previous models.

One of the key parameters of our model is the contact rate for hosts, which clearly depends on the structure of the host population and, provided that the hosts populate a certain region and are moving within it, on the size of that region and the hosts' speed and character of movement. One way of specifying contact rates is to assume that hosts occupy sites in a contact network, and a number of publications have discussed spatially structured host-parasite dynamics basing on that standard assumption. Approaches used include the simulation of probabilistic cellular automata (see e.g. Jeltsch et al. 1997) , the derivation of modified versions of SIS/SIRS differential equations (Lion and Boots 2010) and so-called correlation dynamics (see e.g. van Baalen 2002 and references therein).

However, using contact networks does not seem appropriate in the situations we envisage, where branching process models may be applicable, as the dynamics assumes host populations that are constantly mixing. Moreover, instead of stipulating contact rates (by means of contact networks), for such populations, the natural question that arises is how the contact rate will be determined by more basic parameters such as hosts' speed, the character of their movement etc. In Sect. 5, assuming chaotic independent movements of the hosts, we present a few remarks on how the change in the size of the enclosure a given population of hosts inhabits can affect the "encounter rate" for the hosts -the key parameter of the model describing the "effective density" in the host population.

Finally, in Sect. 6 we consider a multi-stage modification of our model that can be used to analyse farm or aquaculture situations in which there are many enclosures or tanks of animals, and an outbreak of an infectious disease occurs in one of them. We assume that the pathogen might at some stage be transmitted to the next enclosure, where the epidemic process starts anew etc., and discuss possible scenarios of the development of the epidemic in the farm. Section 7 contains a few final remarks concerning biological interpretation of our results.

Assume that we have a large population of hosts that can be infected by parasites of one of two types that will be denoted by T 1 and T 2 . The pathogen types can differ in both their lethality and transmission rate. The numbers of infected hosts at time t are represented by the vector

being the time t number of hosts infected with the type T i pathogen (for simplicity, at any given time, any given host can be infected with one type of the pathogen only). We do not keep track of the number of hosts that remain uninfected (susceptible), assuming instead that this number will remain large enough during the time period for which our mathematical model is intended, so that the dynamics of the process {Z(t), t ≥ 0} to be described do not change over time.

A general description of the model is as follows. We assume that each infected individual lives a random time (which will tend to be shorter when one is infected with the "more lethal" of the two pathogen types). During its lifetime, an infected host can encounter susceptible hosts and, with a probability depending on the type of the pathogen it carries, transmit the parasite to them. The rate of such (random) encounters will be specified by a special parameter that we can vary in order to model changes in the density of the host population.

Finally, we also allow the pathogens to mutate, so that when a host originally infected with T 1 encounters a susceptible host, the latter can become infected with T 2 -type parasites (and the other way around). Now we will present a formal mathematical model. First recall that the exponential distribution with rate (or intensity) α > 0 has density of the form

with mean 1/α, and plays a special role in probability theory due to its unique memoryless property that makes the distribution ubiquitous in the theory of Markov processes. Namely, a random variable τ > 0 modelling, say, the time of the first encounter of a given infected host with a healthy one, has this property if, for any s, t > 0, the conditional probability of the event {τ > s + t} given that τ > s coincides with the probability of {τ > t}:

In words, if, at time s we know that there has been no such encounter, then the conditional distribution (given that information) of the residual random time τ − s till the encounter will be the same as the original distribution of τ . It is obvious that if τ has density (1) then

and so the property (2) is clearly satisfied. An equivalent formulation of the property can be given in terms of the distribution's hazard rate r τ (s) that quantifies the probability that, given there has been no encounter up to time s, there will be one "immediately afterwards": in case a random variable τ has a continuous density p τ , the hazard rate is defined by

and, as h ↓ 0,

where, as usual, o(h) denotes a quantity that vanishes faster than h: o(h)/ h → 0.

It is obvious from (3) and (4) that the hazard rate of a distribution is constant if and only if it is exponential (in that case, the hazard rate simply equals the distribution's rate). In applications, one often uses exponentially distributed random variables to model times between successive events of a particular kind and also lifetimes. This is because, on the one hand, such assumptions make sense from the modelling view point (in a large population, meeting a new individual during a time interval [t, t +h], h > 0, can scarcely depend on one's "life history" prior to time t) and, on the other hand, as the resulting models are usually Markovian, so that the powerful machinery of the theory of Markov processes is applicable.

Our basic model assumptions are as follows:

(a) The initial population contains z i hosts infected with parasites of type T i , i = 1, 2. (b) A host infected with type T i parasites lives a random time which is exponentially distributed with parameter α i > 0, i = 1, 2. Clearly, the pathogen with a higher rate α i will be the more lethal one, as the mean lifetime in that case will be lower. (c) Any infected host can encounter susceptible ones. The time till the first encounter of a given infected host (of any type) with a susceptible host is a random variable exponentially distributed with rate λ. It is clear from the above discussion of the properties of exponential distributions that, at any time t, the residual times till encounters of the current infected hosts with susceptible hosts are all exponential with the same rate λ. A similar observation applies to the residual lifetimes; this ensures that the process {Z(t)} will be Markovian: given the current state of the process, its future evolution does not depend on the past one.

In a modification of the model, one can assume that, at time t, the time till the first encounter of a given infected host (of any type) with a susceptible host has a hazard rate λ(t) which can depend on t. This will enable one to model changes in the density of the host population that occur over time (the higher λ, the more often are encounters, which corresponds to higher host density situations). The process will remain Markovian, but will become time-inhomogeneous. (d) At any encounter with susceptible hosts, a T i -infected host meets only one susceptible host (there can be several such encounters during the host's lifetime). At each such instance, the T i -infected host transmits the parasites to the susceptible one with probability β i (so that, with the complement probability 1 − β i , the encounter will have no consequences for the susceptible host). (e) Mutations T 1 ↔ T 2 are possible. A host infected with T 1 -type pathogens will remain such till its end, but, when transmitting pathogens to a susceptible host during an encounter, the newly infected host will carry T 2 -type pathogens with a probability μ 1 . Likewise, μ 2 denotes the probability that a successful transmission of parasites from a T 2 -infected host to a susceptible one resulted in making the latter T 1 -infected. (f) All the above-mentioned random times (lifetimes, times till the first encounter) are independent of each other.

Of course, the above assumptions oversimplify the real biological processes. There are likely to be several or many strains of pathogen, and the probability of mutation from one to another will vary. We suggest that simplifying the system to two strains is likely to retain the same key dynamic features. Further, the distributions of times until events occur are likely to be approximately exponential as argued earlier, and we do not see any reasons why host survival times and encounter rates should depend in any way on each other. Note that "encounters" between hosts are simply occasions where a pathogen can be transferred between hosts. They do not have to involve physical contact. Thus a pathogen transferred by aerosols might be transferred between pigs that are isolated in neighbouring stalls.

The diagram in Fig. 1 illustrates the "physical" meaning of our assumptions. The two horizontal segments represent the mean lifetimes of hosts infected by the pathogens of our two types: the longer segment (of length 1/α 1 ) corresponds to T 1 which we assume less lethal by stipulating that α 2 > α 1 , whereas the shorter one (of length 1/α 2 ) corresponds to T 2 .

When the host population density is relatively low (say, represented by the value λ = λ , as depicted), the lifetime of T 2 -infected hosts will be too short compared to the mean time 1/λ between encounters to give them a good opportunity to encounter susceptibles and hence further propagate in the host population. One may expect that this will result in the eventual prevalence of T 1 pathogens who have better chance of being transmitted as they live longer. However, if the host population density increases (say, to λ = λ , as depicted), then the more lethal type T 2 may have frequent enough encounters which, combined with its higher transmissibility, can lead to its eventual prevalence. As we will see in Sect. 3, the above argument is confirmed by mathematical analysis.

Assumptions (a)-(f) imply that our process {Z(t)} is actually a two-type time homogeneous Markov branching process in continuous time, see e.g. Section V.7 in Athreya and Ney (1972) . That is, {Z(t)} describes the dynamics of a population consisting of individuals (or "particles") of two types, 1 and 2. The transitions of different particles in the process are assumed to be independent. A particle of type i lives for an exponentially random time with rate a i . At the end of its life it disappears. It can either (i) simply disappear (in terms of our modelling assumptions above, this means that a given T i infected host died having never encountered a susceptible host), or (ii) produce one particle of the same type (meaning: there was an encounter, but no transmission occurred); we think of the "newly produced" particle of type i as just the "old" infected host of type T i who keeps living-note that, due to the memoryless property of the exponential distribution, such an identification of a "new" particle with the "old" one is in agreement with our assumptions (a)-(f) (so that the life of one T i -infected host can actually be represented by a succession of several type i particles, for which reason we do not use T i to denote the type of particles in the branching process), or (iii) produce two particles of the same type i (meaning: there was an encounter and successful transmission, but no mutation; one of the "newly produced" particles is actually the original host, the other represents a newly infected-with the same T i -type pathogen-host), or (iv) produce two particles of different types (meaning: there was an encounter and successful transmission and mutation; one of the "newly produced" particles is actually the original host, the other represents a newly infected-with the pathogen of the other type-host).

To see that both sets of assumptions (a)-(f) and (i)-(iv) describe the same dynamics of {Z(t)}, it suffices to note that in both cases we deal with time-homogeneous continuous time jump Markov processes (which follows from the exponentiality assumptions) and then verify that, choosing suitable parameter values for the second model, one can obtain the same transition rates as for the first one.

To do that, we first observe that assumptions (i)-(iv) imply the branching property which means that, for any s ≥ 0, given Z(s) = (z 1 , z 2 ), the future {Z(s +t), t > 0} of the process will follow the same probability laws as that of the sum of z 1 independent copies of Z(t) starting at time 0 with a single particle of type 1 and z 2 independent copies of Z(t) starting at time 0 with a single particle of type 2. It is clear that the first model [specified by (a)-(f)] has the same property.

Moreover, the branching property implies that to completely describe the evolution of the process, it suffices to specify transition probabilities

from the basic initial states e 1 = (1, 0) and e 2 = (0, 1) for arbitrary small time increments h (for more detail, see Chapter V in Athreya and Ney 1972) . It is obvious that the probability of having more than one transition during a small time interval (0, h) will be o(h), so we just need to consider where a single transition can take the process from a basic state e i according to assumptions (a)-(f) and show that the transitions will have the same rates as for a process specified by (i)-(iv) (for a suitable choice of parameter values).

Suppose that, in our branching process, particles of type i have exponentially distributed random lifetimes with rates a i = α i + λ, i = 1, 2. Moreover, at the end of a type i particle's life, it produces a random number of children (possibly of both types) according to the offspring distributions q i ( j 1 , j 2 ) = Pr a particle of type i gives birth to j 1 particles of type 1 and j 2 ones of type 2 given in Table 1 .

Then, due to independence and (5),

1 (0, 0) = Pr initial type 1 particle dies in (0, h), no children produced = Pr initial type 1 particle dies in (0, h) × q 1 (0, 0) These distributions result in transition probabilities as specified in (a)-(f)

which is clearly the same as the probability for a single T 1 -infected host (from the first model) to die during (0, h). Likewise,

1 (2, 0) = Pr initial type 1 particle dies in (0, h), producing 2 children of type 1 = Pr initial type 1 particle dies in (0, h) × q 1 (2, 0)

The corresponding transition in the first model is as follows: a T 1 -infected host did not die during (0, h), but met a susceptible; the pathogen was transmitted, no mutation occurred. Due to independence, the probability of this will be

Virtually the same argument shows that p (h) 1 (1, 1) = Pr initial type 1 particle dies in (0, h), producing one child of each type = Pr initial type 1 particle dies in (0, h) × q 1 (1, 1)

coincides with the probability for a T 1 -infected host to meet a susceptible and transmit the mutated form of the pathogen. Finally, given Z(0) = e 1 , the most likely state for Z(h) is e 1 . In our branching process this occurs when either the original particle lives through the time interval [with probability 1−(α 1 +λ)h+o(h)] or it dies prior to h producing a single type 1 particle [of which the probability is 

In the first model, to get to this state, the initial T 1 -infected host survives for h time units and either has no encounters with susceptibles or has one, but without transmission of a pathogen. The probability of this is

the same value as for the branching process model. Of course, we could also work out the probability for (1, 0) just as the complement probability

but the presented argument demonstrates the difference between the interpretation of the elements of the two models (T i -infected hosts in the first model and type i particles in the second one are not the same, as we noted earlier). The same calculations are applicable in the case of the initial state e 2 , which leads us to the transition probabilities presented in Table 2 (for h → 0, the additive terms o(h) being omitted for brevity).

is the expected number of type j particles present in the process at time t given that the process started at time 0 with a single particle of type i. Clearly,

the identity matrix, and, using the branching and Markov properties of {Z(t)}, it is easy to demonstrate that {M(t), t ≥ 0} possesses the operator semigroup property:

Indeed, given that Z(0) = (z 1 , z 2 ), the value of Z(s) is just the sum of z 1 independent copies of Z(s) starting at time 0 with a single particle of type 1 and z 2 independent copies of Z(s) starting at time 0 with a single particle of type 2. As the process is time-homogeneous, we infer that

From here, using the Markov property and the double expectation law for conditional expectations, we have

which is equivalent to (6). Relation (6) implies (cf. p. 202 in Athreya and Ney 1972 ) that one has the matrix exponential representation

where A 0 = I and

is the so-called infinitesimal generator of (8) is clearly a solution, as seen from its term-wise differentiation.

Evaluating matrix exponentials is rather straightforward: it basically reduces to calculating the values of the function on the matrix' spectrum (for more detail on functions of matrices, see e.g. Chapter V in Gantmakher 1989) . If, say, A is diagonalisable, so that there exists an invertible matrix Q (with inverse Q −1 : Q −1 Q = QQ −1 = I ) such that

for some σ ± (which is the case in our situation, as we will see below), then clearly A = QDQ −1 ,

and so on:

Thus, to derive the dynamics of the means, we just have to compute the generator A, which can easily be done using Table 2 . Indeed, we infer from the table that, as h → 0,

and similarly for M 2i (h). From here and (9) we immediately obtain

Now solving the characteristic equation det(A−σ I ) = 0 for σ we find the eigenvalues σ ± of A given by σ ± = 1 2 γ 1 + γ 2 ± , = (γ 1 − γ 2 ) 2 + 4δ 1 δ 2 (cf. similar calculations of the threshold parameter for a somewhat different two-type model in Section 8.4 of Mode and Sleeman 2000) , with the respective (right) eigenvectors

The eigenvalues σ ± are clearly different from each other (in any case, this is guaranteed by the fact that σ + is the Perron-Frobenius root for the quasi-positive matrix A, cf. Section A.8 in Thieme 2003) , which ensures that A is diagonalizable and we can take the transformation matrix Q from (10) to be given by (8) and (11) imply that

4 The asymptotic composition of the pathogen types

As σ + > σ − , it is clear from (12) that σ + is the so-called Malthusian parameter of the branching process that determines the long-term behaviour of the branching process means. Note that this "long-term behaviour" still corresponds to the initial stage of the epidemic that we are aiming to model (even more so since the means of the branching population "settle" in that asymptotic composition exponentially fast, as shown below). As we said in Sect. 1, the most interesting for us is the case of supercritical processes for which σ + > 0, implying unbounded exponential growth of the population (unless it becomes extinct at a pretty early stage). Otherwise, the process {Z(t)} would be doomed to die out very soon, so that no epidemic would arise.

It is clear that, in the supercritical case, the ratio of the time t expected number of type 2 particles to that of type 1 particles will be given by

if the process starts with a single type 1 particle, and by

provided that it started with a type 2 particle (recall that u − < 0 and σ + −σ − = > 0). Therefore, using (7), we see that, regardless of the initial state (z 1 , z 2 ), the eventual ratio of the mean number of T 2 -infected hosts to that for T 1 -infected hosts will be one and the same quantity

which is a well-known fact from the theory of multitype branching processes (it follows, for instance, from Theorem 1 on p. 185 in Athreya and Ney 1972 , see also p. 203). The convergence rate is clearly exponential: the remainder term decays as e − t . Note also the obvious facts that R 1 (0) = 0, R 2 (0) = ∞ and that R 1 (t) (R 2 (t)) is an increasing (decreasing) function of t (so that always R 1 (t) < R 2 (t)).

Thus the single value R = R(α 1 , α 2 , β 1 , β 2 , μ 1 , μ 2 , λ) completely specifies the eventual balance of the mean numbers of individuals of different types in our supercritical process (whatever the initial values). This reflects a much deeper result on the long-term behaviour of {Z(t), t ≥ 0}-namely, the fact that, with probability one, the scaled vector e −σ + t Z(t) will converge, as t → ∞, to a non-trivial random vector whose distribution is concentrated on the ray {r v, r ≥ 0} collinear to the (positive) left eigenvector v of A corresponding to the Perron-Frobenius eigenvalue σ + (see e.g. Theorem 2 on p. 206 in Athreya and Ney 1972 and references therein) . This implies that convergence to R holds not only for the ratio of the means, but for the random variables Z 2 (t)/Z 1 (t) as well: if we denote by A the event {Z(t) = 0 for all t > 0}, then

(up to an event of probability zero). In words, this means that either the branching process becomes extinct in finite time or the sizes of the subpopulations of individuals of the two types grow unboundedly in such a way that their ratio tends to R.

Observe also that the above shows how fast the composition of the population will change if the encounter rate λ switches to another value. Suppose that the initial value is λ . As we saw, after some (exponentially short) time, the balance of types in the process will establish around the value R(λ ) ≡ R(α 1 , α 2 , β 1 , β 2 , μ 1 , μ 2 , λ ). Now if the value of λ quickly changes to λ , then the population will re-establish balance at a new level R(λ )-again exponentially fast, with the rate characterized by the new value of (provided, of course, that the process will still be supercritical, i.e. σ + > 0 for the new value of λ).

As we discussed earlier, the increase in the encounter rate λ ought to be beneficial for the parasite type with higher lethality as that increases its chances to spread in the host population. Indeed, we have

since α 2 > α 1 by assumption and it is obvious that |γ 2 − γ 1 | < . Figure 2 displays the dependence of R on λ varying in (0,20), for different levels of the lethality α 2 (left pane) and mutation rate μ 2 (on the right) of the type 2 pathogen. In this example, our model corresponds to many previous models in that transmissibility is greater for the more virulent pathogen. The model shows the expected increase in the frequency of more lethal parasites in a host population when the density of the hosts increases. Observe that the threshold value R = 1 (after which type 2 pathogen dominates the population) may play no critical role: our model is a crude approximation for the initial stage of an epidemic only, so R/(R + 1) will just give the proportion of the carriers of type 2 parasite in the population of infected hosts at the end of that stage.

It is important to note that it is the "splitting" of the single "virulence" characteristic of the parasite into two (lethality and transmissibility) that made such a capability possible: if, say, there is no difference in lethality (α 1 = α 2 ) then, as a simple algebraic Fig. 2 The plots of R as a function of λ ∈ (0, 20). For the fixed common values α 1 = 0.5, μ 1 = 0.2, β 1 = 0.3 and β 2 = 0.6, the left pane displays the plots of R for four different lethality levels α 2 = 1, 1.5, 2 and 2.5 for fixed μ 2 = 0.2 (the lower the value of α 2 , the higher the curve), whereas the right one shows the plots for different mutation rates μ 2 = 0.2, 0.3, 0.4 and 0.6 (the higher the mutation rate, the lower the curve), for fixed α 2 = 2 Fig. 3 Convergence of R i (t) to R as t → ∞. For a common set of parameter values, the plots display the behaviour of R i (t), i = 1, 2, for λ = 2, 6, 10 and 14 (from left to right, top to bottom), with the respective values R ≈ 0.210, 1.076, 1.500 and 1.699. In all the cases, the process is supercritical (σ + > 0) calculation shows, the value of R does not depend on density λ. The last observation we could have actually made earlier, as it follows from the model construction. Figure 3 illustrates the stated exponentially fast convergence of the ratios R i (t) to a common limit R in four situations that have common parameter values α 1 = 0.5, α 2 = 1.5, μ 1 = μ 2 = 0.2, β 1 = 0.3, β 2 = 0.6, but different encounter rates. The plot in the top left corner displays the graphs of R 1 (t) < R 2 (t) in the case when the encounter rate λ = 2 is small enough to allow type 1 parasite to dominate (R ≈ 0.210). On the top right plot, we see that, for λ = 6, there establishes a rough balance (R ≈ 1.076), whereas on the plots in the second raw we see type 2 parasites to gain dominance pretty fast (which is due to higher values of = σ + − σ − ), with the limiting values R ≈ 1.500 and 1.699, respectively. The character of the dependence of R on the transmission probabilities is illustrated in Fig. 4 . For four different values of the encounter rate (λ = 2, 6, 10 and 14), the figure shows the plots of R as a function of (β 1 , β 2 ), restricted to the regions where the process is supercritical (i.e. σ + > 0). The values of the other parameters are α 1 = 0.5, α 2 = 1.5, μ 1 = μ 2 = 0.2. As one could expect, the value of R strongly depends on (β 1 , β 2 ) and is an increasing function of β 2 and a decreasing one of β 1 .

In all four cases presented in Fig. 4 the threshold value R = 1 is exceeded (that is, the more virulent pathogen becomes dominant) only when the transmissibility of type 2 pathogen is greater than that for type 1 (β 2 > β 1 ), so it may appear that inequality is a necessary condition for type 2 to prevail. This, however, is not true: it turns out that a higher mutation rate from type 1 to type 2 can compensate for some lack of transmissibility. Figure 5 shows the plots of R as a function of the mutation probabilities (μ 1 , μ 2 ) ∈ (0, 1) 2 for different encounter rates, all other parameters being fixed and common, with the transmission probability for type 1 being double that for type 2 (β 1 = 2β 2 = 0.4). On the left plot corresponding to λ = 4, the maximum value of R barely exceeds 0.5, whereas on the right one, due to the increase in the encounter rate to λ = 7, not only the supercriticality region is much bigger, but also the maximum value is R = 3. We see type 2 parasite's domination in the region where the mutation rate μ 1 (from type 1) is high enough, while μ 2 is relatively small.

Finally, we turn to the dependence of R on the lethalities α j . As one can easily see, This is quite natural, as the increase in a pathogen type's lethality does not improve its chances to prevail when all other parameters of the model remain unchanged. The character of the dependence is illustrated in Fig. 6 showing R as a function of (α 1 , α 2 ) ∈ (0, 7) 2 , for λ = 2, 6, 10 and 14. Note how the encounter rate λ influences the size of the region where the process is supercritical (thus, for λ = 2 it shrinks to a narrow strip in the (α 1 , α 2 )-domain, corresponding to small values of α 1 ). This is consistent with previous models of virulence evolution-it appears that they all would show that a more virulent pathogen strain would not be competitive if all its other parameters were equal to other strains. The "trade-off" relation between virulence and transmission rate in those models is necessary for more virulent pathogens to dominate. The recent model in Day and Gandon (2006) which also does not assume a tradeoff, produces a similar result.

We do not assume a trade-off between transmission and virulence (i.e. a function linking these characteristics), or that there is a genetic covariance between them. As noted in Day and Gandon (2006) , virulence and transmissibility are affected by both the pathogen and its host. But our model suggests that more virulent strains can dominate only if higher transmissibility or mutation rate (or a combination thereof) favours them. It is this subset of more virulent pathogens that can dominate more easily at high host density.

As we said, our main motivation was to model the effect that a change in "effective density" represented by the parameter λ (the rate at which infected hosts encounter susceptibles) can have on the ratio of the number of hosts infected, say, with T 2 pathogen to that for T 1 -infected hosts. But how does the value of λ relate to physically measurable parameters of the modelled situation -for instance, the density of fish in an aquaculture tank (given the tank size and all other parameters of the model are fixed) or the size of the tank (given the number of fish in it is fixed)? How will λ change if one, for example, "squashes" the same host population to a "world" whose linear dimensions are twice as small as for the original one?

The answer to this type of questions will in general depend on what one assumes about the character of the hosts' movements (and of course, on the pathogen transmission mechanism-but we will not address this aspect in our simple analysis in the section). One of the most popular models for "wandering particles" is the famous Brownian motion process {W (t), t ≥ 0} (see e.g. p. 169 in Karlin and Taylor 1981) , which can be thought of as a continuous analog of a simple (symmetric) random walk. Recall that the Brownian motion is defined as a continuous time process with continuous trajectories that starts at zero at time t = 0 and has independent Gaussian increments: N (0, h) for t, h ≥ 0. One of the key properties of the process is its self-similarity: for any a > 0, one has

i.e. these two processes have the same distribution.

Since the total number of hosts is assumed to be large enough, encounter rates are mostly determined by the "local" characteristics (the density of the population and the dimensionality of the space) of the enclosure and will have little dependence on the "shape of the world". Therefore, for analysis purposes, we will assume in this section that hosts perform independent Brownian motions in a "simple world" S in one, two or three dimensions (starting at some "individual" initial points) and show how the encounter rates change when one "contracts" the original world without changing its shape, i.e. switches from S to the set εS = {εx : x ∈ S}. In this context, the dimensionality can actually be thought of as a crude description of the shape of the world.

Suppose that there are N susceptibles in the population and that the movements of all the hosts are independent of each other. The value of the parameter λ = λ N gives the average number of encounters of a given infected host with susceptible ones per time unit (we assume that N is large enough so that the "conversion" of some susceptibles into infected hosts during our modelling "time horizon" does not significantly change the number N and hence the encounter rate λ). It is clear that λ N = N λ 1 and, moreover, that λ 1 = 1/t * , where t * is the mean time to encounter of our infected host with a given susceptible host. Thus the answer to the question on how the encounter rate will increase if the host density living in a "fixed world" increases is simple: it is just proportional to the number of hosts in the world. However, to understand the effect of the world size change (when we switch from S to εS) on the encounter rate λ = λ(ε), we will have to analyze that effect on the mean time t * = t * (ε).

The one-dimensional case First we assume that S = [0, 1] ⊂ R. Our hosts move according to independent Brownian motions, reflecting from the boundaries 0 and 1 of the set S. This can be formalised by introducing the function

and letting the location of our infected host to be given by H 0 = H 0 (t) = ϕ(H 0 (0) + W 0 (t)) and that of a given susceptible one by H 1 = H 1 (t) = ϕ(H 1 (0)+W 1 (t)), where H i (0) ∈ S are some fixed initial positions and W i are independent standard Brownian motions, i = 0, 1. We say that the hosts have an encounter at time t if H 0 (t) = H 1 (t). Now consider the space εS, where our hosts move according to {ε H i (t), t ≥ 0} and denote by t * (ε) the mean time to encounter of the hosts in this new world, ε > 0. It is obvious from the self-similarity property (13) that t * (ε) = ε 2 t * (1) and therefore

Thus, if our hosts are confined to a one-dimensional world where they wander at random, the encounter rate displays inverse quadratic dependence on the world size: say, halving the "living space" will increase the encounter rate fourfold.

The two-dimensional case To avoid dealing with any boundaries, we assume that our hosts live on the two-dimensional sphere S 2 = {(x, y, z) ∈ R 3 : x 2 + y 2 + z 2 = 1} and that our hosts wander on it at random according to independent spherical Brownian motions {H i (t), t ≥ 0} (see e.g. Karlin and Taylor 1981, Chapter 15, Section 13I) , starting at some fixed distinct points H i (0), i = 0, 1.

In this case, Pr(H 0 (t) = H 1 (t), t ≥ 0) = 1, so we need to modify our definition of encounter. Fix a small enough δ > 0 and define the encounter time of the hosts H i as inf{t > 0 : r (H 0 (t), H 1 (t)) = δ}, where r (·, ·) is the geodesic distance on S 2 . Denote by t * δ = t * δ (1) the mean value of this time (suppressing the dependence of the initial locations H i (0); to simplify the exposition, we deliberately make it somewhat sketchy).

Next we consider the "contracted world" εS 2 , where our hosts wander according to {ε H i (t), t ≥ 0}, but the definition of encounter remains unchanged (the hosts should find themselves within distance δ of each other); the mean time to encounter for this case is denoted by t * δ (ε). Again using self-similarity, we can easily conclude that

However, we want to relate t * δ (ε) to t * δ (1), so it remains to clarify the relationship between t * δ/ε (1) and t * δ (1). It is obvious that t * η = t * η (1) is a decreasing function of η > 0. As we are interested in situations where δ/ε is small (despite the small size of the "world", encounters are still relatively rare), it suffices to find the asymptotic behaviour of t * η as η → 0. To do that, we first observe that analyzing the dynamics of the position of H 0 (t) relative to H 1 (t) shows that finding the mean time when the two points are first within distance η of each other is equivalent to finding the mean time a Brownian particle H * (t) (with an initial position at a distance r (H 0 (0), H 1 (0)) from the "North Pole" P = (0, 0, 1) ∈ S 2 and local diffusion coefficient √ 2 times that for the original spherical Brownian processes) will need to get within distance η of P.

Denoting by V (t) the projection of H * (t) on the z-axis, one can easily see from Itô's formula that {V (t), t ≥ 0} is a diffusion process with the state space [−1, 1] governed by the stochastic differential equation

where W 0 is a standard univariate Brownian motion process and the drift and diffusion coefficients are given by

respectively (see e.g. p. 194 in Matthews 1988 ; for convenience we assumed that H * (t) follows a standard spherical Brownian motion on S 2 which will have no adverse implications for the validity of our analysis). The geodesic distance from H * to P is equal to η iff its projection on the z-axis equals r := cos η, so that we need to find

This can easily be done using the standard technique of the method of differential equations (see e.g. Problem B on p. 192 in Karlin and Taylor 1981) : the function v(z) is the bounded solution to the equation

with the boundary condition v(r ) = 0. Setting u(z) := (z 2 − 1)/4, we notice that u (z) = z/2 and so (15) is equivalent to

Therefore the general solution to (15) is given by

which is bounded on (−1, r ) iff c 1 = 1. Now using the boundary condition at z = r to find c 2 leads to v(z) = 4 ln

To find the asymptotic behaviour of v(z) as η → 0, we use cos η = 1−η 2 /2(1+o(1)) to obtain that, for a fixed initial value z,

This means that, for fixed initial positions of H 0 and H 1 , we have t * η = (c + o(1))| ln η| as η → 0. Together with (14) this yields, for small δ/ε,

That is, the encounter rate behaves as

In the two-dimensional case, it might be more natural to relate the encounter rate not to the linear dimensions of the enclosure, but rather to its area (proportional to ε 2 ). The above formula shows that, in this case, the encounter rate in a "shrinking" world still grows somewhat faster than the inverse proportion of its area, the latter giving the density of hosts.

The three-dimensional case

We still have (14), and an analysis similar to the one carried out in the two-dimensional case shows that now t * η = (c + o(1))η −1 as η → 0 (cf. p. 195 in Matthews 1988) , so that t * δ (ε) = ε 3 t * δ (1) (1 + o(1) ). That is,

We see that, in the three-dimensional case (assuming that the hosts wander according to three-dimensional Brownian motions), the encounter rate in a "shrinking" world is inversely proportional to its volume. That is, in this case the rate is proportional to the density of hosts. To summarise the above analysis, we observe that the encounter rate λ = λ(ε) grows rather fast when the linear dimensions (specified by the parameter ε) of the hosts' "world" diminish. In the three-dimensional case, λ is inversely proportional to the volume per host, in the two dimensional case it grows slightly faster than the reciprocal of the area per host, while in the one-dimensional case the growth rate of λ is inversely proportional to the square of the size of a host's share of the enclosure. This indicates that not only effective density per se, but also the shape of the enclosure can be an important factor leading to an epidemic. Thus the nature of the enclosures in which animals are kept can be an important factor in determining the progress and nature of an epidemic.

In this section we will consider an aggregate model for situations in which there are several populations of hosts that exist in originally isolated enclosures. The situation we model here could be viewed as a simple instance of a meta-population model (we refer the reader to Hanski and Gaggiotti 2004 as an extensive source on metapopulation modelling). But in our scenario there is no feedback from later enclosures to earlier ones, and we do not consider any dynamics of the hosts. Thus we find it more useful to consider it as a chain of epidemic episodes.

In intensive animal husbandry, there are often many pens, paddocks or sheds enclosing animals at a farm or many fish tanks at an aquaculture facility. We consider the scenario where initially, one of the enclosures becomes infected with a single type pathogen. This can give rise to a "local epidemic" in the infected enclosure, which can be modelled using our processes from Sect. 2 (assuming that we have supercriticality: σ + > 0). The original pathogen may also mutate to become more or less lethal.

We will initially assume that it may mutate to a more lethal type 2 pathogen (α 2 = r α for some r > 1, α 1 = α). That new pathogen type can also have different transmissibility and mutation rate, but, to make our model as simple as possible, we will assume for the time being that it differs from type 1 pathogen in lethality only, all other parameters being unchanged. Denoting them simply by β and μ, we see that the Malthusian parameter for that process is given by

At some point, the infection might be transmitted to the next enclosure (say, by a worker in a farm situation). We assume that this occurs after the epidemic has gone through the initial stage, and so the ratio of the numbers of hosts infected with different pathogen type can be assumed equal to ρ 0 , where ρ k ≡ R(r k α, r k+1 α, β, β, μ, μ, λ) . The transmitted pathogen is chosen at random, so that the probability of transmitting the one with lethality α (denote this event by A) is 1/(1 + ρ 0 ), while the one with lethality r α is transmitted with probability ρ 0 /(1 + ρ 0 ). This, in turn, may lead to a local epidemic in the new enclosure: we again assume the possibility of mutation to a more lethal pathogen (so that now we will have α 1 = α, α 2 = r α if the event A occurred, and α 1 = r α, α 2 = r 2 α otherwise), and to have an epidemic we again need σ + = σ + (α 1 , α 2 ) > 0 (now for the new set of parameters α 1 , α 2 ). Once the epidemic has established itself in the second enclosure (and the balance of pathogen types has stabilized around the respective R-value), the next step is the transmission of the disease (by means of a random mechanism of the same type as in the first instance) to the next enclosure, and so on.

Scenarios of this type have been encountered often where once a disease is recognized in a herd, animals in the infected enclosure are removed or killed, but the disease is subsequently found in other herds, for example the spread of foot and mouth disease among herds in Taiwan, which was related to herd size and the number of herds in a province (Gerbier 1999) . Of course, biosecurity measures are intended to prevent such transmission between enclosures, but often the need for diligence is learned after the event.

It is easily seen from (16) that

and also that ∂ ∂α R(α, r α) < 0. Thus, if the lethality of the pathogen will keep increasing, the Malthusian parameter of the model will eventually drop below zero, and then the epidemic will collapse. More specifically, setting k * ≡ inf{k ≥ 0 : σ + (r k α, r k+1 α) < 0}, we see that

It is clear that the transition of the disease from enclosure to enclosure according to the above scheme is described by a discrete time Markov chain {X n }, the "time" n having the meaning of the number of steps (i.e. enclosures infected), X n representing the level of lethality of the pathogens in the nth infected enclosure: we set X n = k if (α 1 , α 2 ) = (r k α, r k+1 α) in the enclosure. Thus the state space of the Markov chain is {0, 1, 2, . . .} and the only nonzero entries in the transition matrix P = [p j,k ] j,k≥0 of the chain are

Further, in view of (18), one can assume that once the Markov chain {X n } has reached the state k * , the epidemic becomes unsustainable, and hence there will be no further transmission of the disease to other enclosures. So we can truncate our state space to {0, 1, 2, . . . , k * }, which results in a finite decomposable Markov chain with a single absorbing state k * . Whatever the current state of the chain, at the next step it can either stay at it or move to the right, the transition probabilities forming the matrix

where T is the k * ×k * substochastic matrix formed by the first k * rows and k * columns of P , r = (0, . . . , 0, p k * −1,k * ) ∈ R k * + and denotes transposition. The (random) number of steps T the chain will need to reach the absorbing state is nothing else but the total number of enclosures that will be affected by the epidemic prior to its collapse. Using our model, we can easily find the distribution of T .

Indeed, using the standard approach to solving such problems (see e.g. p. 80 in Kijima 1997), we note that as the state k * is absorbing, we have Pr(T ≤ n| X 0 = 0) = q (n) 0,k * , where q (n) jk are the n-step transition probabilities:

[q (n) jk ] = Q n = T n r n 0 1 , r n = (I + T + · · · + T n−1 )r , so that for the probability mass vector function

Of course, we are only interested in the first entry of the vector f (n).

To compute the mean and higher moments of T one can use the generating function

In particular, since d dz f * (z) = r I − zT −1 + rzT I − zT −2 , we find that

where the last equality follows from the obvious observation that r I − T −1 = f * (1) = 1 ≡ (1, . . . , 1) ∈ R k * + . A possible objection to the above simple aggregate model is that pathogens will not always mutate to become more lethal. The model can be further generalized by allowing, within each enclosure, mutations of our pathogen not only in the direction of higher lethality, but also in the opposite direction. So we will first have to generalize our basic model from Sect. 2 to a three-type branching process, assuming that, if an enclosure is infected with a pathogen with lethality α, then the pathogen can mutate to ones with lethalities r −1 α and r α, respectively, where, as before, r > 1 is a fixed number (all other parameters being assumed equal for the different types of pathogens). Mathematically, analyzing such processes is essentially equivalent to what we did in Sects. 2-4, although all the closed-form expressions will become much more complicated, and so we will not give much technical detail for brevity's sake. The main assertions concerning the asymptotic behaviour of the branching process will remain true: there will exist a limiting balance of types in the supercritical case (denote the shares of the different pathogens by π j = π j (α), j = 1, 2, 3, j π j = 1), which can be found from the generator A of the semigroup of the mean matrices, and the almost sure convergence of the process scaled by e −σ + t (as before, σ + denotes the Perron-Frobenius root of A) to a limiting random vector will hold.

In the multi-stage model, we start with initial infection of one of the enclosures with a pathogen with lethality α. That leads to an epidemic (provided, of course, that σ + > 0) in which pathogens of three types will be present, with lethalities given by the vector (α 1 , α 2 , α 3 ) = (r −1 α, α, r α). The next enclosure to be infected will receive a pathogen chosen at random from those present in the first infected enclosure, and so it will have lethality α j with probability π j , j = 1, 2, 3, and so on.

Observe that the triplets of lethalities (α 1 , α 2 , α 3 ) characterizing the pathogens present in a given enclosure in our system will all be of the form (x, r x, r 2 x) for some x > 0, i.e. lying on a common ray L with the direction vector (1, r, r 2 ). Therefore we will again have a basically "univariate" Markov chain {X n } showing what pathogens can be present in different enclosures, X n = k meaning that the nth affected enclosure was initially infected with the pathogen of lethality r k α, k ∈ Z (and in this case there can also be pathogens with lethalities r k−1 α and r k+1 α in that enclosure), assuming that X 0 = 0 (as α is the lethality of the pathogen that was initially introduced into our system). The original state space for the process will be Z ≡ {. . . , −1, 0, 1, . . .}, which is infinite in both directions. At each step, the value of the chain can remain unchanged (the interpretation being that the pathogen transmitted to the next enclosure had the same lethality as the one with which the epidemic started in the current one) or can either decrease or increase by one (that is, the transmitted pathogen would have lethality values equal to r −1 or r times the current one, respectively).

Further, as before, one can show that ∂ ∂ x σ + (x, r x, r 2 x) < 0, so that, moving along the ray L in the "positive direction", we will eventually enter the subcriticality region for the branching process, where the basic reproductive number will be less than one. Therefore, at this stage the epidemic will collapse, and hence we again can "truncate" the state space for {X n } -now to {. . . , −1, 0, 1, . . . , k * }, where k * is an absorbing state that has the same meaning as above.

The variety of behaviours that such a chain can display will be somewhat wider than for our first multi-stage model. The dynamics of the chain are determined by the behaviour of the mean step values

In particular, the parameters of the model can be such that the above quantities will be negative. Then absorbtion at k * occurs with probability less than one, while on the complement event the chain will drift away in the negative direction, which corresponds to the disease "fading", when the pathogen's lethality vanishes, and so on.

We conclude with remarks concerning possible biological interpretation of our results. First, our model results are generally consistent with previous models of virulence, but we have focused on the initial stages of epidemics in an animal husbandry context, and used a different approach and assumptions. The mathematical models we have presented show that, at the beginning stages of any epidemics that arise in such situations, the density of hosts is an important factor in determining whether more or less lethal strains of the pathogen will predominate. Many authors have noted that selection for virulence differs during the initial (spreading) stage of an epidemic (Frank 1996; Day and Proulx 2004; Bull and Ebert 2008; André and Hochberg 2005) , and others have noted that high host densities are likely to lead to increased virulence (Ewald 1994; Gandon et al. 2001; Bull and Ebert 2008) .

Many previous models (Ewald 1994; Frank 1996; Day 2003) have assumed that more rapid production of copies of the pathogen inside the host would both increase the likelihood of transmission to new hosts and also shorten the life of the host. There has been considerable criticism of this "tradeoff hypothesis", although there are limited data supporting it (de Roode et al. 2008; Fraser et al. 2007; Bolker et al. 2010 ), but it has now been realised (Bull and Ebert 2008) that it is simply not relevant to the initial stages of an epidemic. Nevertheless, our model shows that increased lethality alone will not lead to dominance by a pathogen strain-although more lethal strains can invade more easily at high densities and they will dominate more rapidly, provided they are favoured by some factor such as higher transmission or mutation rates.

Our models show that the ratio of infections by more and less lethal pathogens stabilizes very fast, so that even if measures to prevent further spread of the disease are put in place as soon as an outbreak is identified (such as elimination of the animals in an enclosure), the relative frequency of pathogen types is likely to have changed before action is taken. Thus a high density of animals on farms or mariculture facilities will rapidly lead to the dominance of more lethal strains of pathogens if these can enter the farm and mutations occur to produce more lethal variants. An example of this process has recently been described in the mariculture of fish in Pulkinnen et al. (2010) . The problem has been recognized in intensive poultry production (Slingenbergh and Glibert 2010) , and the identification of more virulent trains of Marek's disease in chickens (Witter 1997 ) may also be an example of this process. It also seems possible that the advent of very virulent strains of Spanish Flu in August 1918 in US military camps, troop ships, and European disembarkation camps, some time after the start of the epidemic may be linked to high densities of soldiers in these camps and troop ships and the difficulty in implementing quarantine in these overcrowded camps (Mathews et al. 2009 ), although there are other plausible explanations, such as overcrowded conditions leading to higher infection doses.

Thus the increasing density at which animals are kept (intensification of the industry) should be considered as a risk factor for the evolution of more virulent diseases. Various authors have pointed to this issue (Biggs 1985; Gandon et al. 2001; André and Hochberg 2005) but in many cases attention has been focused on the possibility that vaccination-especially where the vaccine reduces pathogen replication-selects for higher virulence. We contend that it is likely that the higher transmission of pathogens at high densities that allows more virulent pathogens (which may replicate faster) to become dominant in spite of such vaccines.

Assuming a chaotic pattern of movement of animals inside the enclosure where they are kept (as modelled by independent Brownian motion processes), we have modelled how the effective stock density affects their encounter rate (the key parameter of our model) and thus influences what pathogen type will predominate. The effect increases more rapidly with density where an enclosure is most appropriately measured as area (cattle feedlots, paddocks) than if the appropriate measure is a volume (fish tanks). We assumed constant mixing so that the rate of encounter with susceptible hosts is not materially changed in the first stages of the epidemic. These assumptions will limit the situations in which this result may apply to cases where animals mix freely and rapidly in the enclosure, but our result may also apply if there is a large area around an infected individual where others may become infected by a pathogen, because it is dispersed in air, water or over time [as a long-lived propagule (Biggs 1985) that can disperse or survive until restocking (Bull and Ebert 2008) ].

In agricultural systems, humans have been maintaining animals and plants at very high densities for about 10,000 years, and humans themselves began living at high local densities in villages, towns and cities from about the same time (Crawford 2007) . Densities of both humans and domesticated organisms have probably increased continually since then, but the densities at which some domestic animals are kept appear to have increased very rapidly in the last few decades (see e.g. Fraser 2005) , and the recent rapid development of large scale commercial aquaculture involves the maintenance of very high densities of a whole new set of marine species that will have their own pathogens (see Harvell et al. 2004 Harvell et al. , 1999 Bergh 2007) .

The ratio of pathogen types will clearly also depend on the transmissibility of pathogen strains. We have focused on differences in lethality between pathogen strains, because pathogen strains with increased transmissibility would obviously become more prevalent. It is interesting to note that our model showed that, even if the trans-missibility of pathogen of one type is lower than that for the other, the former pathogen can still prevail provided it is favoured by mutation. Increased density of animal hosts may itself increase transmissibility, due to stress on the hosts caused by crowding (Sniezko 1974) . Increased transmissibility might in turn lead to the evolution of increased lethality (Ewald 1994) .

Modern animal husbandry often involves a large number of separate enclosures, each containing a large number of animals at very high densities. Once a disease is detected in an enclosure, farmers would either use antibiotics, or sacrifice or remove the animals, but pathogens may be carried between enclosures by various mechanisms, depending on the type of pathogen and the biosecurity practices followed. A multistage version of our model for this situation suggests that if pathogens are transferred a number of times, then the evolution of more lethal pathogens may be very rapid, but the increase in lethality will eventually lead to the epidemic becoming unsustainable (hosts dying too fast to be able to transit the pathogen).

We suggest that the outcomes predicted by the mathematical models discussed in the present paper can carry important messages for animal husbandry, where there are strong commercial incentives to increase the densities of animals in enclosures to very high levels, and often very large numbers of enclosures are built in a single farm.

Risk factors for the evolutionary emergence of pathogens

Virulence evolution and the trade-off hypothesis: history, current state of affairs and the future

An introduction to stochastic epidemic models

Coevolution of hosts and parasites

Infectious diseases of humans: dynamics and control

Stochastic epidemic models and their statistical analysis

New York Ball FG, Donelly P (1995) Strong approximations for epidemic models

Some evolutionary stochastic processes

An introduction to stochastic processes

The effect of heterogeneity on the spread of disease

The dual myths of the healthy wild fish and the unhealthy farmed fish

Essai d'une nouvelle analyse de la mortalité causée par la petit vérole et des advanteges de l'inoculation pour la prévenir

Comparison of populations whose growth can be described by a branching stochastic process

Infectious animal disease and its control

Transient virulence of emerging pathogens

Invasion thresholds and the evolution of non-equilibrium virulence

On the evolution of virulence and the relationship between various measures of mortality

Virulence evolution and the timing of disease life-history events

Insights from Price's equation into evolutionary epidemiology

A general theory for the evolutionary dynamics of virlence

Virulence-transmission trade-offs and population divergence in virulence in a naturally occurring butterfly parasite

Mathematical models for infectious disease statistics

The evolution and expression of parasite virulence

Host-parasite relations, vectors, and the evolution of disease severity

Evolution of infectious disease

Animal welfare and the intensification of animal production: an alternative interpretation. UN Food and Agriculture Organization

Variation in HIV-1 set-point viral load: epidemiological analysis and an evolutionary hypothesis

Imperfect vaccines and the evolution of pathogen virulence

Effect of animal density on FMD spread. European Commission for the Control of Footand-Mouth Disease Report

A simple model of epidemics with pathogen mutation

On the statistical measure of infectiousness

Isolation and characterization of viruses related to the SARS coronavirus from animals in southern China

Branching processes: variation, growth, and extinction of populations

Ecology, genetics and evolution of metapopulations

Emerging marine diseases-climate links and anthropogenic factors

The rising tide of ocean diseases: unsolved problems and research priorities

Extinction times in multitype Markov branching processes

A thousand and one epidemic models

Branching processes with biological applications

Pattern formation triggered by rare events: lessons from the spread of rabies

Deterministic and stochastic epidemics in closed popultaions

A contribution to the mathematcial theory of epidemics

Markov processes for stochastic modeling. Chapman and Hall, London Knolle H (1989) Host density and the evolution of parasite virulence

The evolution of virulence in parasites and pathogens: reconciliation between two competing hypotheses

Are parasites 'prudent' in space?

Virulence and transmissibility of pathogens: what is the relationship

Qualitative behavior of stochastic epidemics

Transmission rates and the evolution of pathogenicity

Understanding influenza transmission, immunity and pandemic threats

Covering problems for Brownian motion on spheres

Applications of mathematics to medical problems

Modeling and real-time prediction of classical swine fever epidemics

Intensive fish farming and the evolution of pathogen virulence: the case of coumnaris disease in Finland

Do old and new forms of poultry go together

The effects of environmental stress on outbreaks of infectious diseases of fishes

Contact networks and the evolution of virulence

Increased virulence of Marek's disease virus field isolates

The outcome of a stochastic epidemic-a note on Bailey's paper

How do pathogen evolution and host heterogeneity interact in disease emergence?

Yersin A (1894) La peste bubonique á Hong Kong

Acknowledgments K. Borovkov was supported by the ARC Centre of Excellence for Mathematics and Statistics of Complex Systems (MASCOS) and ARC Discovery Grant DP120102398. T. Rice was partially supported by MASCOS. The original motivation for this study arose from conversations of R. Day with Jeremy Prince and Harry Gorfine. The authors are grateful to the referee whose thoughtful comments helped to improve the presentation of the paper.