key: cord-0004771-xfi3p601
authors: Trapman, Pieter; Meester, Ronald; Heesterbeek, Hans
title: A branching model for the spread of infectious animal diseases in varying environments
date: 2004-03-03
journal: J Math Biol
DOI: 10.1007/s00285-004-0267-5
sha: 6353f50ebc75a4f88524d38c0311250f692e3b7e
doc_id: 4771
cord_uid: xfi3p601

This paper is concerned with a stochastic model, describing outbreaks of infectious diseases that have potentially great animal or human health consequences, and which can result in such severe economic losses that immediate sets of measures need to be taken to curb the spread. During an outbreak of such a disease, the environment that the infectious agent experiences is therefore changing due to the subsequent control measures taken. In our model, we introduce a general branching process in a changing (but not random) environment. With this branching process, we estimate the probability of extinction and the expected number of infected individuals for different control measures. We also use this branching process to calculate the generating function of the number of infected individuals at any given moment. The model and methods are designed using important infections of farmed animals, such as classical swine fever, foot-and-mouth disease and avian influenza as motivating examples, but have a wider application, for example to emerging human infections that lead to strict quarantine of cases and suspected cases (e.g. SARS) and contact and movement restrictions.

Recent outbreaks of infectious diseases of animals (e.g. classical swine fever (CSF), foot and mouth disease (FMD) and avian influenza (AI)) in Western Europe have had great impact on the economy, public life and animal health and welfare in the countries involved. During such an outbreak one would like to be able to compare the effectiveness of proposed control measures in, for example, their ability to reduce the expected final size and the expected duration of the outbreak. Typical for the strategies aimed at stopping outbreaks of important diseases of farm animals, is that infected herds are removed from the population by culling upon detection. A second characteristic is that due to increasing quantity and quality of the imposed control measures, the environment that the infectious agent experiences, is changing. By this we mean that consecutive measures can make, for example, contact opportunities between herds different in different phases of the outbreak, or can make the infectious period, or rate with which infectivity is produced, differ for farms infected at different times.

In most cases, mathematical methods for computing outbreak characteristics such as expected final size and expected duration, assume a constant environment in that the control measures are not compounded in time and do not lead to changes in the rates that govern epidemic spread (see e.g. [2] ). In this paper we aim to develop stochastic methods, based on branching processes, which allow us to compare the effectiveness of control strategies during such outbreaks in situations where the environment is varying because of changes in subsequent measures of control.

Much work has already been done to describe the spread of classical swine fever (see e.g. [10, 11, 1] ) and foot and mouth disease (see e.g. [5, 6, 9, 8] ). In this paper, we will model the spread of infections in a much more analytic way than is done in earlier models [5, 6, 9, 1] . We use an iterative method that computes properties of the spread, like the probability of a major outbreak of the infection and the final size of an epidemic (i.e. the total number of infected herds) very efficiently. Furthermore we can derive some properties of the duration of the epidemic. We allow for different types of herd.

In our model, it is essential that once the infection in a herd is detected, the whole herd will be culled. Therefore, our main interest is the number of infected herds, but the number of infected animals in an infective herd is important for determining the infectivity of a herd and the distribution of the detection times. We therefore model the spread of the infection on two levels, namely the spread of the infections within a herd and the spread of infection between herds; both are described by a stochastic process.

We use a special branching process to describe the spread of the infection among herds. The parameters of this branching process depend on the time since the infection of the herd and on the environment, which is determined by the real time. Using branching processes to describe epidemics is of course not new (see e.g. [4] ), but no theory exists that gives short-term predictions (as opposed to asymptotics) for general branching processes in varying environments, with an age-dependent birth rate.

For computations, it is necessary that after a certain moment the environment is constant. To achieve this we assume that after some time no new measures will be taken and the effects of all measures taken in the past will either be constant or absent. In other words, although the values of the parameters may differ from those before certain measures were taken, the values are assumed to be constant after a given moment in time.

Our model can be used to predict the effects of various control measures and strategies during an ongoing outbreak. Meester et al. gave a method to estimate the parameters from the data available during an outbreak, [10] . We consider measures like single vaccination of all herds of a certain type, a total transport ban or killing of animals just after birth. Some of these measures cause a varying environment. The fraction susceptible animals in a herd or the fraction susceptible herds of the total number of herds may be varying and so the infection rate may vary in time as well. In fighting outbreaks, additional measures will be implemented as soon as present measures turn out to be insufficient. A change in measure will change the values of the parameters. We assume that the changes in the environment are deterministic.

We develop the theory using a classical swine fever outbreak in herds of pigs as motivating example throughout. In our paper we use the same input data as Klinkenberg et al. [1] .

As mentioned in the introduction, we first need to model the spread of the infection within one herd, since the infectivity of an infective herd, and also the time at which the infection is detected in a certain herd, depend on the number of infected animals in that herd. We use t for the time elapsed since the first measures were implemented. The variable τ is used for describing the spread of the infection within the herds, it is the time since infection of a particular herd and therefore relative. The τ -clock starts ticking at the moment the first animal in the herd becomes infected. From now on, we will call τ the "infection-age" or "age" of the herd.

We use only four types of disease-related parameters:

µ: The recovery rate of individual animals in a herd. λ: The infection rate of individual animals within a herd. α: The per capita detection rate of infected animals within a herd. β: The rate at which one infected animal infects susceptible herds.

We assume that as soon as an infection at a particular herd is detected, the whole herd will be culled instantaneously. Further we assume that the rate of detection is proportional to the number of infected animals at infection-age τ , I (τ ), i.e. the detection rate is αI (τ ), where α does not depend on the infection-age. We assume that the infection in a herd develops as an autonomous process until detection. The number of infected animals in a herd therefore depends only on the infection-age τ , and not on the absolute time t.

We describe the number of infective animals by an ordinary birth and death process, writing p i (τ ) for the probability of i infective animals in a herd at infection-age τ and δ ij for the Kronecker Delta function. We have, (see e.g. [3] ):

Solving these differential equations leads to:

where r = λ − µ and R = λ µ is the reproduction ratio. We assume λ > µ, otherwise the infection will typically only cause a minor outbreak within a herd, and will therefore typically not infect other herds.

Hence, conditioned on the event that the epidemic in a herd does not go extinct before age τ , i.e. I (τ ) > 0, I (τ ) is geometrically distributed with parameter (1 − Rp 0 (τ )) = R−1 Re rτ −1 , which is small for large τ . A geometric random variable with small parameter can be approximated by an exponential random variable with the same parameter. Therefore, for large τ , I (τ ) can be approximated byH (τ ), whereH (τ ) is an exponential random variable with

where H is an exponential random variable with parameter R−1 R . For large τ , the term 1 R is negligible compared to e rτ . So we use the approximation

where the equality is in the distributional sense.

We can interpret this approximation as follows. The random variable H represents the random character of the start of the outbreak in a herd, when only a few animals are infective. If the disease does not go extinct, then after the initial phase, there are many infected animals in the herd, each of which causes an independent number of new infections per time unit. According to the law of large numbers, this means that the growth rate is eventually almost deterministic. We also use this approximation for small τ .

We use the word infective for animals or herds that are able to spread the infection. We use the notation I (τ ; h) for the number of infective animals in a particular herd, with H = h given.

1. We assume that the infection and recovery rate of individual infected animals are independent of time and age. 2. It is possible to use alternatives for the birth and death process to describe the spread within a herd, for example the contact process. The contact process is appropriate when we have a situation where all animals are positioned in a row and do not change position. Each animal can only infect its two nearest neighbours. We assume that recovered animals with two infective neighbours are re-infected immediately. Therefore, we only consider the animals at the edge of a row of infective animals. We still assume that the recovery rate is µ. An infective animal infects each of its susceptible neighbours with rate λ. We have that

Let the number of infectives be I (τ ). Then E(I (τ )|I (τ ) > 0) = 2rτ + o(τ ). Furthermore, we can show that Var(I (τ )|I (τ ) > 0) = o(τ 2 ). This implies that for I (τ ) > 0 we have

So, in contrast to the birth and death model, we may approximate the random variable I (τ ) for large τ by a deterministic variable 2rτ . This makes computations much easier in this case. 3 . In the original model, the number of animals in one herd is assumed to be very large compared to the number of infected animals. Therefore, we assume that the contact rate between infected and susceptible animals is constant, and hence so is the birth rate. From data of outbreaks of CSF in the past, we can see that the number of infected animals until detection of the infection within the herd, is small compared to the total number of animals in the herd, so these assumptions seem justified in this case [11] . 4. The within-herd infection and recovery parameters λ and µ can be measured experimentally or from data of past or on-going outbreaks. The parameter α depends on the development of symptoms of infected animals and how attentive the farmers are. Therefore, in reality this α will change at the first detection of the infection in the country (or in neighbouring countries), due to higher awareness of farmers and veterinarians. It is very difficult to estimate α for the period before the time of the first detection. For the time after the first detection, we can estimate α from data of past outbreaks in the same area or try to estimate this parameter during the ongoing epidemic. Using data from past outbreaks is dangerous, because the characteristics of the virus and of the farming practice may have changed. Estimation of the parameters, during an on-going outbreak is done by Meester et al. [10] . This method has some problems, e.g. the time necessary to get enough data for a reliable estimate. Another problem is that in [10] the infectivity of a herd does not depend on the "age" of the herd. This independence of age is essential for the estimations made. For our model it is not necessary that α is constant in time. It is possible to extend the model and use α(t) instead of α. 5. We assume that culling is the only measure which influences the within-herd spread of the infection. Vaccination is assumed to show no effect on the spread in the herd. For CSF this assumption is justified by the fact that vaccination will lead to immunity only after two weeks. Therefore, during this first two weeks the spread of the infection is not affected by this measure. For the time after these two weeks, we assume that we can use the same speed of propagation of the infection, for computational reasons. The approximation leads to results that are conservative, that is, too pessimistic. 6. We do not take characteristics of individual animals into account that might cause the individuals to differ in infectivity, susceptibility or contact pattern.

Often age and type of species can have a substantial influence. For CSF this implies we do not distinguish between the age of animals in one herd. The detection rate and infectivity of herds with many young animals does not significantly differ from the detection rate and infectivity of herds with especially older animals [11, 1] . 7. We also use the approximation I (τ ) = H e rτ for small τ . This is not correct, but due to the small number of infective animals at small τ , the probability that the disease is detected at small τ is small too. We will see later that the number of infections in that period is as well, so the overall influence of the events while the herd is 'young' is probably not so large.

In this section, we consider classical swine fever as a concrete example. We distinguish between two types of farms: multipliers (m) are,roughly speaking, farms where young piglets are born, and finishers (f ) are farms that buy piglets and fatten them. p m denotes the fraction multipliers of the total number of herds and p f is the fraction finishers. We assume that within either of these types, a birth and death model (with the same parameters for both types) describes the within-herd spread. The infectivity per non-transport contact of both types of herds develops in the same way too. However, transport contacts are only allowed from multipliers to finishers. Therefore, we take a larger infectivity rate for contacts from multipliers to finishers.

We define A ξ (t, τ ; h) as the infectivity of a herd at time t, while the herd was infected τ time units ago and with H = h, where ξ is a two-dimensional vector, denoting the two types of herds involved in the contact, so that ξ can be ff , f m, mm or mf . When no measures are implemented, A ξ (t, τ ; h), with h a given realisation of H , is proportional to the number of infected animals in the herd.

Here β ξ is a constant depending only on ξ . Because the non-transport contacts all happen at the same rate we can define β m := β f m = β mm . Further we write β f := β ff and β mf is β f plus some additional term for infections caused by transport of piglets from multipliers to finishers. Because all non-transport infections happen at the same rate the ratio β m /β f is exactly the ratio of the multipliers to the finishers. Note that the infectivity does not depend on the absolute time t.

We define β by β := β m + β f , hence β m = p m β and β f = p f β. In order to consider transport contacts we also define β mf = p f β + β tr , where β tr is the proportionality factor of the part of the infection rate that is due to transport contacts. There is no term p f in front of β tr , because all transport contacts are from multipliers to finishers. We assume that the total number of herds is very large; in our computations, we assume it to be infinite.

We already know that the detection-rate is given by αI (τ ; h) = αhe rτ . From this we can deduce, for given h, the probability p nd (τ ; h) that an infected herd of age τ is not yet detected is

In the case where no measures are implemented, the expected number of infected multipliers by one infective multiplier up to age τ , for given h, µ mm (τ ; h), is given by:

To see this, we first consider only one type of herd. We prove the following proposition:

Proof. The infection and the detection rate are proportional to the number of infective animals in a herd. The ratio of infection rate and detection rate is given by β/α. We call the detection of the herd and the infections by this herd events. The probability that an event is an infection is β α+β and a detection α α+β . Therefore, the number of events, including detection, is described by a geometric random variable with parameter α α+β . Hence the direct offspring distribution does not depend on the size at the start of the process. Moreover, the same holds for the offspring of this direct offspring.

In the same way we can prove that the size of the future offspring of an infective herd is independent of its age.

To prove (3), consider two types of herds. We note that 1 added to the number of multipliers infected by one infective multiplier is given by a geometric random variable with parameter α α+β m . (1 is added because the final event will be the only detection, and the number of events is described by a geometric random variable.) So the expected number of future infections of an infective multiplier is β m α , at all times. The probability that an infective herd is not yet detected at age τ is given by e − α r h(e rτ −1) . From this and Proposition 1 we deduce that the expected number of infections after age τ is given by β m α e − α r h(e rτ −1) . By subtracting the expected number of infections after age τ from the expected total number of infections by one herd we get the expected number of infections until age τ .

In the same way we can deduce µ mf =

). Now consider the probability p f kl that a finisher infects k multipliers and l finishers. All events (infections of multipliers, infections of finishers and detection of an infected herd) happen at a rate proportional to the number of infective animals in the infective herd. Therefore, the rates are all proportional to each other. First, we only consider infections and detections, and we do not yet consider the different types of herds infected. Detection occurs with rate αhe rτ and infection occurs with rate β m he rτ + β f he rτ = βhe rτ . As in the proof of Proposition 1 we can describe the total number of events,D say, by an ordinary geometric random variable with parameter α α+β . So the probability that n + 1 events occur, i.e. n herds are infected by one finisher, is

If in total n herds are infected by one finisher, we know by the lack of memory property that the number of infected multipliers, N f m , is binomially distributed with parameters n and β m β . That is,

We can deduce

in the same way. Note that we did not need the distribution of the random variable H , we only needed the proportionality factors of the parameters.

In the special case where infection and detection rates are proportional and these rates are known for the time before the first detection, we can also give the distribution function of the number of infective herds at the time of the first detection. This distribution function is the same as the distribution function of the number of direct infections by one herd, because it still holds that an event is a detection with probability α α+β and the probabilities of infections of finishers and multipliers are respectively β f α+β and β m α+β . We are also interested in the probability that a herd infects k multipliers and l finishers, before the herd reaches age τ . We consider a finisher. We write P(N f m (τ ; h) = k, N ff (τ ; h) = l) for this probability when the detection age τ d and H = h are given. The infections before age τ d occur independently of each other and have the lack of memory property. Therefore, for τ ≤ τ d the times of infections of the different types of herds are described by an inhomogeneous Poisson-process with rate β ξ he rτ . So the number of infections until age τ is Poisson distributed with parameter hβ ξ τ 0 e rs ds = hβ ξ (e rτ − 1). Now we see:

Note that now we know the life length distribution, the expected number of infections by a herd up to age τ and the underlying Galton-Watson process of the branching process, we can use the general theory of branching processes to determine some other properties, like the expected duration of the outbreak. ( [4] , chapters 2 and 6).

1. We use a branching idea to describe the spread of the infection among different herds. For this approach, we assume that one herd has contacts with many other herds. In our model we do not take spatial distribution of the herds into account. As long as there are many herds in an area, the local exhaustion of susceptible herds can be ignored. If the outbreak is in an advanced state, however, local depletion of susceptible herds (either by the epidemic progression or by so-called pre-emptive culling) will certainly play a role. As long as there are relatively few infected herds in a neighbourhood, (i.e. a group of herds that have contacts with each other) we assume that the infection rate does not directly depend on this number of infected herds. Measures like ring-culling (culling of all herds within a certain distance of an infected herd) and ring vaccination cannot be considered in our model. 2. We distinguish between two types of contacts, transport contacts and indirect contacts. We have assumed that the indirect contacts are the same for all herds. Transport contacts are only from multipliers to finishers. Indirect contacts include transport of the virus by wind, by visits from an infective herd to a susceptible herd etc. Transport contacts are transports of infected animals from a multiplier to a finisher. 3. We simplify the model by assuming that the measures are implemented at the time of the introduction of the virus in the population (or, in other words, that the infection is detected immediately). It would be desirable to implement the measures from the moment of first detection. It is difficult, however, to estimate parameters like α and β for the time before the first detection, when awareness is not yet heightened by announcement of the outbreak, and when increased hygienic measures on farms are not yet taken. If we know the value of the parameters of the model for the period between introduction and first detection, we are able to estimate the probability of extinction and the expected final size of such a situation. However, these calculations still require that we use the information from our simplified model, where measures are started directly upon introduction. 4. The function µ mm (τ ; h) is not required for our calculations of the final size and the probability of extinction. The reason why we include this function in our paper is that this function is interesting in its own right. It gives the expected offspring τ time units after an infective herd is infected itself. In practical applications, it is possible that by contact tracing some herd is suspected of being infected τ time-units ago. This µ mm (τ ; h) is then the expected number of infected multipliers by the suspected herd, if that suspected herd is a multiplier and if it is really infected. 5. The proportionality of the infection rate to the detection rate is essential in this section. Due to this property, we can find an underlying ordinary Galton-Watson process. Without this proportionality, we may ask whether it is possible to give a nice expression for the generating function. We also lose the independence of h in the generating function, which gives some computational problems. 6. Knowing the generating function of the number of infected herds at the first detection is important because if in some way it is possible to estimate parameters for the time before the first detection, we can use the distribution function of the number of infective herds at the time measures are implemented. We do not have the 'ages' of the infective herds at the first detection, but we can find a worst-case-scenario, by finding what age of a herd at t = 0 leads to the biggest offspring. If the detection and infection rate does not depend in the same way on the number of infected animals, we cannot deduce the distribution function at the time of the first detection in this way.

In the previous section, we considered the model for the spread of an infectious disease in non-varying environments. We heavily used the proportionality of the infection and detection rate of herds. This proportionality does not hold in a varying environment. Therefore, we need a different approach. We consider only one type of herd; the multi-type model is a straightforward generalisation of this. We assume that the spread within a herd is not influenced by the state of the environment. The detection rate only depends on the number of infected animals in a herd and is written as αI (τ ; h), where I (τ ; h) is the number of infected animals as described in Section 2. In this section, we assume that the within-herd spread is described by a birth-and-death process. So I (τ ; h) = he rτ . We can easily deal with other descriptions of the within-herd spread.

The infection rate may change due to control measures. Some measures lead to an infection rate that only changes finitely many times, other measures cause a continuously varying infection rate. We assume that in either case the environment is constant after some given time, t = T . That is the moment that measures have no more added value and therefore do not lead to new changes in the values of the model parameters. The infection rate is given by

where φ(t) describes the effects of the measures implemented. We assume that φ(t) is deterministic. Because the environment does not influence the within-herd spread, φ(t) does not depend on the age τ . For the rest of this section it is assumed that φ(t) = 1 for t ≥ T .

The assumption of constant environment after some given time is often realistic. With a vaccination for instance, we know the moment that all vaccinated animals no longer live. In other control measures, like a transport ban, we can vary the time T and compute the effects.

We are looking for the probability of extinction of the infection, the expected final size and a generating function for the number of infected herds at a certain moment. To do this we use a discrete approximation of φ(t), so we have only a finite number of changes. Because the environment is non-varying after a certain moment, we can use the ordinary theory of branching processes to get all interesting properties for herds infected after t = T . By using backward iteration we will find the relevant properties of the epidemic for a herd infected in another interval, because these properties only depend on what happens in the intervals after the interval of infection. We can compute these properties for the herd infected at t = 0.

Because an infected herd will almost surely be detected in finite time, the probability of extinction of the 'progeny' of an infected herd x, is equal to the probability of extinction of the progeny of all the herds infected by x. The probability of extinction of the "progeny" of a herd infected at time t, is denoted by q(t). So:

Here, the times t i are the times that herds are infected by the herd infected at time t; the expectation is over these random times of infection. The empty product is defined as 1. In order to make computations possible we use a discrete time approximation. We divide the positive real line into N + 1 intervals, labeled 1, 2, . . . , N + 1. The

The final interval N + 1 is (T , ∞); in this final interval all the model parameters are constant. The time t = 0, the moment of the first infection, is not included in one of the intervals, we treat this point in our notation as interval 0.

If the function φ(t) is discontinuous we may choose another discretisation, so that the discontinuities are on the boundaries of intervals. It is not essential that all intervals have the same length.

For t in interval i, 1 ≤ i ≤ N, q(t) and φ(t) are approximately constant. In our "discrete time model" we will write q(i) and φ(i) for the value of q and φ in interval i. For instance, we can take φ(i) to be the value of φ(t) at the midpoint of interval i. Because we accept discontinuities in φ(t), this function may differ significantly for neighbouring intervals. Now, with n(i, j ), the number of infections in interval j , due to one herd infected in interval i, we can write (4) as:

where the expectation is over the numbers of infections in each interval. q(0) is the probability of extinction of the progeny of the herd infected at time t = 0. We assume that all infections and detections take place on the midpoint of the interval, except of course for the events in the final interval, because there we can compute everything explicitly. Given that a herd is not yet detected and H = h, the probability of infection in a certain interval by that herd does not depend on the number of infections in previous intervals. Consider a herd infected in interval i, denote by D(i; h) = k the event that this particular herd is detected in interval k, given H = h. After detection no further infections occur. Define n(i, l, k; h) as the number of infections in interval l due to a particular herd infected in interval i, while H = h and the interval of detection k, are given. Now, by using independence we have (writing P H for the distribution function of H )

Here the expectation inside the integral depends on h, contrary to the non-varying environment case. Note that the number of infections in a certain interval is not independent of the interval of detection. It does not make any difference for the number of infections whether a herd is detected shortly after the considered interval or very long after that, but it does make a difference whether the detection is in the considered interval itself. So p(n; i, l, k; h) := P(n(i, l, k; h) = n) depends on the detection interval k. For computational reasons we suppose in this discrete model that for i ≤ N, p(0; i, i, k; h) = 1. We write p det (i, k; h) for P(D(i; h) = k). We are interested in q := q(0). We can easily compute all these probabilities in our model.

Because we assume a birth-and-death process for the within-herd spread, we have to take the random character of H into account. Doing this leads to the following formulae:

· · · E(q(k) n(i,k,k;h) )p det (i, k, h)dh.

Here we have conditioned on the time of detection and on h. We can rewrite this formula as:

Here q(i) depends on q(l) for l > i. As mentioned before, we can compute q(N + 1) by the ordinary theory of branching processes. We use backward iteration to compute q(0).

Note that we also know the probability of extinction for herds infected after time t = 0. If we estimate (for example by contact tracing) the moment of infection of a certain herd, we can use this information to improve predictions.

In almost the same way we can calculate the expected final size of the epidemic i.e. the expected total number G of infected herds. This will also be done for only one type of herd.

If q < 1, there is a positive probability for the final size to be infinite and therefore the expected final size will always be infinite; so to have the possibility of a finite expected final size we assume q = 1. We denote by G(t) the size of the progeny of a particular herd, infected at time t (including this ancestor) and write G(0) = G. The expected number of herds in the progeny of herd X, including X itself, is 1 plus the sum of the expected size of the progenies of all herds infected by this X, i.e.

Where the t i 's are again the random times at which infections by the element, infected on time t, occur. The empty sum is defined to be equal to 0.

In the same way as we deduced the formulae for q(i) we deduce the formulae for E(G(i)) in the discrete approximation. We writeḠ(i) for E(G(i)):

For α > β,Ḡ(N + 1) is known from the ordinary theory of branching processes (see e.g. [4] ): for φ(N + 1) = 1,Ḡ(N + 1) = α α−β .

We will give the generating function for the number of infective and infected herds at any given moment. Here we will consider the generating function for the infective and infected herds at time t = T , the time after which the parameter values are considered to be constant. Remember that the age of an infective herd only influences the number of infective animals in that herd. From Proposition 1 in Section 3.1, we know that the expected direct offspring of a herd, born after t = T , is independent of the age of that herd, given that the herd is not yet detected at that time. Furthermore we know that the size of the offspring (infected after time t = T ) of a herd, infective at that time, does not depend on the offspring of other herds that are infective at time t = T . So for the distribution of the number of infections after time t = T , only the distribution of the number of infective herds at that time is important.

By using Proposition 1 and the theory for ordinary branching processes, we can compute everything we want to know. We will determine the distribution of X, the number of infective herds at time t = T . X i is the number of infective herds at time T , from the progeny of one particular herd infected in interval i. (Here the herd itself is a part of its own progeny.) We denote X 0 by X. We also use Y i for the number of infected herds in the progeny of a particular herd infected in interval i, that is detected before time t = T . The generating function of the distribution of X i and Y i is

We writeg i (s 1 , s 2 ; h) = E(s X i 1 s Y i 2 |h). We again assume that a herd does not infect other herds in the interval wherein it becomes infected itself. So: For interval N we have:

Now we can determineg 0 (s 1 , s 2 ) point wise. Note that we only need to computẽ g 0 (s 1 , s 2 ) in m + 1 points to give a good approximation of the first m derivatives of this function for a certain point. With this derivatives we are able to compute the first m moments of the size at time t = T or to approximate P(X = n) for all n ≤ m.

With this generating function we can also compute the probability of extinction and the expected final size of the epidemic. We can use the same reasoning as before. The progenies of all herds, infective at time t = T have to go extinct. This occurs with probability q X 0 . So the probability of extinction is ∞ k=0 P(X 0 = k)(q(N + 1)) k =g 0 (q(N + 1), 1). Note that with these s 1 and s 2 the model is exactly the same as the model in Section 4.2. For the expected final size we add the expected number of infected herds, already detected, to the expected size of the progeny of all the infective herds at time t = T (again including the herds infective at that time). This is d

. We can only use this property if after time t = T the infection and detection rate are proportional to each other, because otherwise we need to know the ages of the infective herds at time t = T .

In a non-varying environment, we know for q = 1 the "speed" at which P(Z(t) > 0) decreases, for t → ∞, where Z(t) is the number of infective individuals at time t, descending from one individual infected at time t = 0 [4] . So we can give an upper bound for the probability the disease is already extinct at a certain time. We do this by assuming that all infective herds at time t = T are infected at that time. So all herd at time T have age τ = 0. Now with the notation S t for the probability that the progeny of a herd infected at time T , will not survive until time T + t and with Z(t) for the number of infected herds at time t, we have that:

1. By using this iterative method, we can compute the expected final size of an outbreak very fast, but for computing any higher moments of this final size, we need substantially more computational effort. By using simulations, it is possible to estimate these higher moments too. (see [1] ) 2. The lack of memory property is very important for our computations. We used it to write the formula of q(i) in a 'convenient' form, with expectations in front of every q(l) n . The speed of computations also heavily depends on the lack of memory property of the infection rate. In each interval the number of infections by one herd is Poisson distributed. We can easily integrate h out of the formula for q(i), for Poisson distributed numbers of infections. (With the given distribution of H ).

In this section we use the Dutch classical swine fever epidemic of 1997 as example.

We use the same data as Klinkenberg et al. in [1] . However, Klinkenberg et al. used simulations and for every simulation the parameters were (pseudo-)randomly chosen from the distributions of the parameters, while we used only estimations of parameters for our computations. We consider the following set of control measures:

A Total transport prohibition, B Killing of young piglets, in combination with a breeding ban, C Vaccination of all piglets (not sows) at multiplier herds, followed by recurrent vaccination of newborn piglets, D Single vaccination of all pigs at finishing herds, E Vaccination of piglets on arrival at finishing herds.

Choosing a (combination of) measure(s) is called a control strategy or scenario. The effects of different scenarios on the fraction of infective animals in a multiplier, in a finisher and the possibility of transport-infections (φ m (t), φ f (t) and φ tr (t) respectively) are given in Table 1 (This table is adapted from [1] ). For example, the 0 for φ tr in strategy A means that infection is not possible by transport contacts. As another example, φ f (t) = t/100 for t ≤ 100 in strategy D means that at time t a fraction t/100 of the animals in a finisher is infective. Note that the first 12 strategies are in a constant environment. Therefore, for those strategies we can compute properties of the spread of the disease directly. The probability of extinction is given in Table 2 .

We also compute the expected final size of an epidemic; see Table 3 . We only need to compute this for the strategies with almost sure extinction. The expected number of infected multipliers, while the initially infected herd was a finisher is denoted by G f m . In a similar way we define G mm , G mf , G ff . The last column of Table 3 is the expected number of infected herds, when initially five multipliers and five finishers were infected. Note that we cannot compare the expected final size with the results of Klinkenberg et al. (Table 6 in [1] ), because they use the median of 1000 simulations, not the mean. The results of Klinkenberg et al. are given in Table 4 . The 95% confidence intervals of the final size are also taken from [1] and were estimated by simulation.

We used the generating function for the number of infected and infective herds at time t = T , the point in time after which the parameter values are assumed to be constant, to estimate the size at time T (which varies for different strategies). Table 1 . The effect of different strategies on φ f (t), φ m (t) and φ tr (t) 1 . We estimate the first two moments of the size at time T for the four pure strategies in a varying environment B, C, D, E. (Table 5 ) Especially the results of B and C are of interest because they will give an idea of the variance of the final size of a herd.

Next we give covariance matrices (for different strategies and different initially infected herds) of the number of infected multipliers not yet detected, the number of infected multipliers already detected, the number of infected finishers not yet detected and the number of finishers already detected, respectively. We denote by V ar i (J ) the covariance matrix for an initially infected herd of type i, while the strategy is J . By comparing the values in these matrices to the results in Table 5 , we know how much trust we may put in the expected sizes. Up to now we used exact parameters. In reality we cannot estimate the parameters exactly. In order to get some insight into how the computed properties depend on the values of the different parameters, we varied one parameter while keeping the other parameters constant. The results for the strategies B and D are given in Figures 1 and 2 . Note that the scales on the axes differ for the different varying parameters. On the x-axes we put the ratio between the value of the parameter and the point estimator of that parameter. We see that for these strategies, the parameters r and R have relatively little influence on the computed quantities, while α and β do have significant influence.

In a constant environment the only thing that matters is the ratio β α . For varying environments we varied α, while keeping β α and β tr α constant. The results are varying, but not very much.

1. We cannot estimate the parameters exactly. Because the expected final size, the probability of extinction and the generating function for the number of infective herds at a certain moment heavily depend on some of the parameters, we cannot really use the computed values in a quantitative way. But we can use them to compare different scenarios, while using the same parameter values. 2. We computed the probability of extinction for the progeny of one infective herd, but in reality we may have more than one infected herd at the moment the first measures are implemented, so the probability of extinction may be much less than the computed q. We assumed that different infected herds infect other herds independently of each other. If we simplify the model by assuming that all herds at time t = 0 have age 0, we can estimate the expected final size by the number of initially infected herds times the estimated final size of one infected herd. In the same way we can estimate the probability of extinction by the computed probability for one herd, to the power of the number of initially infected herds. From previous outbreaks of FMD, CSF and AI in the Netherlands, one can see that as a rule several farms are already infected at the moment the first case is suspected or confirmed. 3. The results are very much the same as the results of Klinkenberg et al. but we do not need simulations to estimate the properties. Also, our method is much faster than using simulations. In Table 6 of [1] it is suggested that the expected final size, with initially five infected multipliers and five infected finishers, is estimated by simulations. In reality the median of 1000 simulations is used in that paper. This explains why our results differ on that point. 4. Note that strategy ABC gives an extremely high expected final size. Klinkenberg et al. [1] did not even give this final size. This is because for the given parameters the process is very near to a process with a positive probability of surviving i.e. a positive probability of an infinite final size. 5. Using the expected size at t = T and the variance of this size, we see that the expected final size is not very informative in some scenarios. Large variance may imply that a large set of "numbers of infected herds" have a not-ignorable probability, even if we know the exact parameters. 6. The computed covariance matrices do have values of different order in them.

Consider for instance strategy B. The variance of the number of finishers and multipliers already detected is much higher than the variance of the number of finishers and multipliers that are still alive at time t = T . This is because the infection rate decreases for this strategy. At the start of the process one infective herd may infect several other herds with a relatively high probability, while after some time infections become rare. This is why the expected number of infective herds at t = T is much less than the expected number of infected herds. 7. For constant environments we are interested in the variance of the final size.

Note that in all of the strategies with almost sure extinction in a constant environment, only one type of herd can be infected, so we only have to deal with one type of herd. We use the underlying Galton-Watson process to deduce the variance of the final size (see Section 2.11 of [4] ). If we take m for the expected number of direct infections by one infective herd, the expected final size is 1 1−m and the variance is m(1+m) (1−m) 2 . Note that for a large expected final size, the variance will be of the square order of the final size. Due to this large variance, we cannot give exact quantitative predictions about the final size of an epidemic. This also indicates that there is a large intrinsic uncertainty in the problem. 8. If the infection rates are increasing in time, the expected final size will decrease for increasing r, while for decreasing infection rates the expected final size will increase for increasing r. This is because r can be seen as a parameter describing the speed of the process in a constant environment. Large values of r correspond to short generation lengths compared to small r. So, if the infection rate is increasing in time, the n-th generation for small values of r will have a larger infection rate than the n-th generation for large r.

The effect of decreasing R, the parameter describing the random effects at the start of the within-herd spread, is less clear from formulae and definitions. From the figures we can only see that R has relative small effect on the expected final size and the probability of extinction. 9. Some parameters heavily influence the expected final size and the probability of extinction. We have already seen that in a constant environment r and R do not influence these quantities. For predictions the ratios β : β tr : α are most important. So the estimation effort is best devoted towards estimating these ratios.

By using a stochastic model, we could estimate the probability of extinction and the expectation of the final size of an epidemic in a varying environment. We only need that the environment is not varying anymore after some given time. We use a branching process in varying environments, depending on the age of the particles. In the constant environment the probability of extinction and the final size are known. By using an iterative process, we computed this probability and size for the time the environment is still varying. By using generating functions, we can compute many important properties, like the moments of the final size, and a lower bound for the probability that the process has gone extinct after a given time. These generating functions are very useful especially if in a constant environment the expected number of infections after a certain age does not depend on that age. This only holds if the infection and detection rate are in direct proportion, for the non-varying environment.

It is difficult to estimate the parameters for the computations. Some of the properties computed in this model, like the probability of extinction heavily depend on the parameters α and β, describing the infection and detection rate. The model presented in this paper may be useful to compare the effect of different measures, but it is very dangerous to use this model for absolute quantitative predictions. We used the model to describe the spread of classical swine fever. It is worth investigating if the same model, with other parameters, may be used to describe the spread of other animal diseases, with culling at detection, or even human diseases, which lead to strict quarantine of detected infected individuals (and suspected cases) and contact and movement restrictions. Emerging infections such as SARS are possible examples of this.

Quantification of the Effects of Control Strategies on Classical Swine Fever Epidemics

Mathematical Epidemiology of Infectious Diseases

Probability and Random Processes

Branching Processes with Biological Applications

The Foot-and-Mouth Epidemic in Great Britain: Pattern of Spread and Impact of Interventions

Transmission Intensity and Impact of Control Policies on the Foot and Mouth Epidemic in Great-Britain

Branching Processes with Deteriorating Random Environments

The Role of Mathematical Modelling in the Control of the 2001 FMD Epidemic in the UK

Dynamics of the 2001 UK Foot and Mouth Epidemic -Dispersal in a Heterogeneous Landscape

Modelling and Prediction of Classical Swine Fever Epidemics

Quantification of the Transmission of Classical Swine Fever Virus between Herds during the 1997-1998 Epidemic in The Netherlands