key: cord-0177175-ridg3vis
authors: Borovkov, Konstantin; Day, Robert; Rice, Timothy
title: High host density favors greater virulence: a model of parasite-host dynamics based on multi-type branching processes
date: 2011-08-31
journal: nan
DOI: nan
sha: ad58d6a70c13eb218d487794e5568d9d82cce157
doc_id: 177175
cord_uid: ridg3vis

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. Instead of using a single virulence characteristic which is typical of most mathematical models for infectious diseases, our model uses a combination of two characteristics: lethality and transmissibility. This makes the model capable of reproducing the empirically observed fact that the increase in the host density can lead to the prevalence of the more virulent pathogen type. 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.

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 D. Bernoulli's paper [13] (see also pp. 2-6 in [16] ). A major step in further development of deterministic modelling was the "threshold theorem" [38] 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 [43] , applying the "law of mass action" suggested in [13] to probabilistic description of the epidemics. Stochastic models became much more popular after Bartlett [8] formulated a model for the general stochastic epidemic by analogy with the deterministic model from [38] . May and Anderson [47] introduced the first model of evolutionary change in pathogen virulence in 1983, 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 [47, 19, 44, 40, 23] , 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 [2, 3, 6, 10, 16, 49] and also to papers [4, 18, 34] .

In natural systems, the density of host organisms declines when a virulent pathogen produces an epidemic that kills host rapidly. The models of virulence have made clear that this leads to selection for less virulent strains, 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. Thus one can 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 [29] or when the plague bacterium, Yersina pestis, has been transferred to humans from rodents [56] .

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 [38] on pp. 11 and 29 in [16] ) insights into the dynamics of pathogen populations where the virulence of the pathogens can change due to mutation. More specifically, we are interested in analyzing the effects that changes in host density may have on the virulence of the pathogen. 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 [15] . 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. [22] , 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 [32, 31, 12] ). Thus the question arises as to whether high host densities are likely to lead to selection for more virulent strains of pathogens.

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 [14] and also refer to [9, 37, 55] . 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 [11, 48, 33] .

Branching processes are relatively simple and well-understood mathematical models; for treatments of the theory of branching processes, see e.g. [35, 5] , a recent addition to the monographic literature on branching processes prepared specifically for biological audience being [30] . 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. [42, 7, 34] ). 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.

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 [49] ).

An additional argument for this approach is that in agriculture or aquaculture situations, when a threat of an epidemic outburst arises, drastic measures are usually taken: the affected part of the population may be isolated or even destroyed. These measures will substantially change the way pathogens can pass between hosts, and hence can make the standard stochastic epidemic models inapplicable beyond the initial stage of the epidemic anyway.

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.

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 is 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. One of the advantages of the branching process model is that one can easily incorporate in it two different characteristics of pathogenes' virulence: their lethality (which can be described by the mean survival time of an infected host) and transmissibility (specifying the probability of an infected host infecting a susceptible one on their contact), instead of representing them together by a single quantity termed simply "virulence" (cf. e.g. [1] ). For recent discussions of the interplay between these two characteristics, see [27, 41] .

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. The underlying simple continuous time Markov model of a two-type branching process is presented in Section 2. The analysis of the model and derivation of the dynamics for the mean functions are given in Section 3 and used in Section 4 to establish the eventual composition of the pathogen population. It turns out that our model is capable of reproducing the above-mentioned phenomenon of shifting the overall parasite's lethality in response to increased density of the hosts. This is actually possible because of using two different parameters to characterise the virulence of the pathogen.

In Section 5 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 Section 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 Figure 1 : Mean lifetimes of infected hosts and mean times to encounter. 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 Section 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 [5] . 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. 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 by the following table: Table 1 . Offspring distributions in the branching process model (i)-(iv).

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)

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

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). 

Consider the matrix of mean values

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 [5] ) that one has the matrix exponential representation

where A 0 = I and A = lim

is the so-called infinitesimal generator of {M (t)}. Indeed, from (6) (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. Chaper V in [24] ). 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 ,

(11) 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 

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 [53] ), which ensures that A is diagonalizable and we can take the transformation matrix Q from (10) to be given by

Hence (8) and (11) imply that

4 The eventual composition of the population 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 process means. As we said in Section 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 [5] , see also p. 203, ibid.). 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 longterm 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 {rv, 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 [5] 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 lim t→∞ Z 2 (t) Z 1 (t) = R on A (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 it gains more time to spread in the host population. Indeed, we have

since α 2 > α 1 by assumption and it is obvious that |γ 2 − γ 1 | < ∆. Fig. 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. This example demonstrates that our model is capable of reproducing the growth of 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 a proportion of the carriers of type 2 parasite in the population of the infected hosts at the end of that stage. Even a relatively law value of that quantity may be fatal for the population.

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 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, resp.

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 the four cases presented in Fig. 4 the threshold value R = 1 is exceeded only when the transmissibility of type 2 pathogen is greater than that for type 1 (β 2 > β 1 ), so it may appear that 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 pathogene 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 ). 5 On the effects of the change in enclosure size on the encounter rate As we said, our main motivation was to model the effect that the 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 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 [36] ), 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 though of as a crude description of the 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

, 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

and that our hosts wander on it at random according to independent spherical Brownian motions {H i (t), t ≥ 0} (see e.g. [36] , 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 [46] ; 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 [36] ): 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, v(z) = 8| ln η| + O(1).

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 proportional to 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 [46] 

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 inverse 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.

There may be many pens with animals at a farm or many fish tanks at an aquaculture facility. 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 Section 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 common. Denoting them simply by β and µ, we see that the Malthusian parameter for that process is given by

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 infection is transmitted to the next enclosure (say, by a worker in a farm situation). 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 [25] . 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 σ + (r k α, r k+1 α) < 0 for all k ≥ k * . (18) In particular, since

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 Section 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 Sections 2-4, although all the closed-from 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, rx, 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, rx, 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 E X n+1 − X n X n = k) = π 3 (r k α) − π 1 (r k α), k < k * .

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 will conclude with a few remarks concerning possible biological interpretation of our results. The mathematical models we presented show that, at the beginning stages of any epidemics that arise in situations where animals (or humans) live in enclosures, the density of hosts is an important factor in determining whether more or less lethal strains of the pathogen will predominate. 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. Other key factors that can also contribute to the evolution of more lethal strains and their spread in the host population include transmissibility and mutation rates. Mutation rates are known to vary greatly between pathogens such as the flu virus, and others such as trematode worms (liver flukes, shistosomes etc., which reproduce more slowly).

Our models show that an increase in the 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 [50] . The problem has been recognized in intensive poultry production [51] . and the identification of more virulent trains of Marburg's disease in chickens [54] may also be an example of this process, and it seems possible that the advent of very virulent strains of bird flu during the Spanish Flu epidemic may be linked to high densities of soldiers in demobilization camps and troop ships, although these men may have been very susceptible due to their poor condition [45] .

Thus the density at which animals are kept should be considered as a risk factor for the evolution of more lethal diseases. Assuming chaotic character of movement of animals inside the enclosure where they are kept (as modelled by independent Brownian motion processes), we discussed how their effective density affects the key parameter of our model specifying the hosts' encounter rate and hence eventually determining what pathogen type will predominate.

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 transmissibility of pathogen of one type is lower than that for the other, the former pathogen can still prevail provided it is favoured by mutation. Density of the animal hosts may itself increase transmissibility, due to stress on the hosts caused by crowding [52] . Increased transmissibility might in turn lead to the evolution of increased lethality [20] . An interesting question is whether increased lethality would be linked in most cases to increased transmissibility. Some previous models [20, 21, 17] 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. We did not make any assumptions of that kind for our model.

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 usually 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.

Population dynamics of microparasites and their invertebrate hosts

Infectious Diseases of Humans: Dynamics and Control

Stochastic Epidemic Models and Their Statistical Analysis

An introduction to stochastic epidemic models

Branching Processes

The Mathematical Theory of Infetious Diseases and its Applications

Strong approximations for epidemic models

Some evolutionary stochastic processes

An Introduction to Stochastic Processes

Analysis of Infectious Disease Data

The effect of heterogeneity on the spread of desease

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

Deadly companions: how microbes shaped our history

Epidemic Modelling: An Introduction. Cambrdige

Virulence evolution and the timing of disease life-history events

Mathematical models for infectious disease statistics

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

Evolution of infectious disease

Models of parasite virulence

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

Imperfect vaccines and the evolution of pathogen virulence

The theory of matrices

Effect of animal density on FMD spread. European Commission for the Control of Foot-and-Mouth Disease Report. Research Group of the Standing Technical Committee

Basic methods for modeling the invasion and spread of contagious diseases

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

Emerging marine diseasesclimate 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

A Second Course in Stochastic Processes

Deterministic and stochastic epidemics in closed popultaions

A contribution to the mathematcial theory of epidemics

Markov Processes for Stochastic Modeling

Host density and the evolution of parasite virulence

Virulence and transmissibility of pathogens: what is the relationship

Qualitative behavior of stochastic epidemics

Applications of mathematics to medical problems

Transmission rates and the evolution of pathogenicity

Understanding influenza transmission, immunity and pandemic threats

Covering problems for Brownian motion on spheres

Epidemiology and Genetics in the Coevolution of Parasites and Hosts

Modeling and real-time prediction of classical swine fever epidemics

Stochastic processes in epidemiology

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

Mathematics in population biology

Increased virulence of Marek's disease virus field isolates

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

(1894) La peste buboniqueá Hong Kong

Acknowledgments. K. Borovkov was supported by the ARC Centre of Excellence for Mathematics and Statistics of Complex Systems.

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 areFurther, 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 matrixwhere 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 [39] ), 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:so that for the probability mass vector function f (n) = {f j (n), j = 0, . . . , k * − 1}, f j (n) = Pr(T = n| X 0 = j), one obtains f (n) = r n − r n−1 = r(T n−1 ) ⊤ , n ≥ 1.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 f * (z) ≡ ∞ n=1 z n f (n) = rz I − zT ⊤ −1 , |z| ≤ 1.