- Email: [email protected]
Stochastic dynamics of immunity in small populations: A general framework
Mathematical Biosciences 200 (2006) 28–43 www.elsevier.com/locate/mbs
Stochastic dynamics of immunity in small populations: A general framework A.-F. Viet
a,b,*
, G.F. Medley
a
a
b
Ecology and Epidemiology Group, Department of Biological Sciences, University of Warwick, Coventry CV4 7AL, United Kingdom Unit of Animal Health Management, Veterinary School – INRA, BP 40706, 44307 Nantes Cedex 03, France Received 9 August 2004; received in revised form 8 December 2005; accepted 22 December 2005 Available online 8 February 2006
Abstract Assessment of immunological status is a powerful tool in the surveillance and control of infectious pathogens in livestock and human populations. The distribution of immunity levels in the population provides information on time and age dependent transmission. A stochastic model is developed for a livestock population which relates the dynamics of the distribution of immunity levels at the population level to those of pathogen transmission. A general model with K immunity level categories is first proposed, taking into account the increase of the immunity level due to an infection or a re-exposure, the decrease of the immunity level with time since infection or exposure, and the effect of immunity level on the susceptibility and the infectivity of individuals. Numerical results are presented in the particular cases with K = 2 and K = 3 immunity level categories. We demonstrate that for a given distribution of the immunity levels at the population level, the model can be used to identify quantities such as most likely periods of time since introduction of infection. We discuss this approach in relation to analysis of serological data. 2006 Elsevier Inc. All rights reserved. Keywords: Epidemic; Immunity; Stochastic model; Population; Transmission dynamics; Seroepidemiology
*
Corresponding author. Address: Unit of Animal Health Management, Veterinary School – INRA, BP 40706, 44307 Nantes Cedex 03, France. Tel.: +33 (0)240 687 806; fax: +33 (0)240 687 768. E-mail address: [email protected] (A.-F. Viet). 0025-5564/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2005.12.027
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
29
1. Introduction The rational control of an infectious disease depends on data to assess the incidence and prevalence. Such data usually consist of the health status at time t of the population based either (1) on the record of disease occurrence at individual level before t, or (2) on the presence of an indicator of a previous infection at individual level at time t. The first approach requires that the disease is recognisable and accurately recorded for each individual. However, most infections are asymptomatic or show unspecific signs, for example bovine-viral diarrhoea virus (BVDV) in cattle [1], and respiratory syncytial virus (RSV) in humans [2]. Consequently, specific diagnosis frequently relies on serology (i.e. finding specific antibodies to the infection within an individual). The second approach relies on surveying the population for pathogen-specific antibody levels as an indicator of a previous infection. Such data give a measure of the extent of previous infection and, indirectly, of the prevalence at the time of the survey. For examples of the use of such serology see Bahri et al. [3] and Maloo et al. [4]. Cross-sectional and longitudinal surveys of pathogen-specific antibody levels provide valuable data for assessment of infection dynamics, especially when combined with information on age to provide age serological profiles (e.g. [5–7]). The immunity level of an individual is determined by measurement of the concentration of specific antibodies to the pathogen in a sample of serum or, more recently, oral fluid (‘saliva’). For most clinical needs, it is sufficient to categorise individuals into two classes: positive (and protected) and negative (and susceptible to infection). As assessment of antibody concentration is quantitative (titre), this requires that the original continuous measurement is reclassified into a binary score based on a predefined or relative cut-off value. However, it has recently been realised that the distribution of titres itself contains valuable information [8,9]. Because titre generally increases during or immediately after an infection or a re-exposure and decreases with time since exposure, the antibody level potentially contains information on the time since infection of the individual. Thus, distribution of specific individual antibody levels contains indirect information on the transmission dynamics of the pathogen. For example, if the pathogen has been present, but is now absent, a few individuals would have a high immunity level and most individuals a low level, whereas, with the pathogen present, most individuals would have a high immunity level. Previous work has considered the relation between pathogen transmission and the distribution of the immunity level at the population level. Some researchers aimed at interpreting such observed distributions in terms of pathogen spread or of time since infection of the population (e.g. [1,10]). In contrast, Glass and Grenfell [11] and Woolhouse et al. [12] introduced antibody levels to consider waning protection in childhood disease and against infection with foot and mouth disease virus, respectively, and the impact of vaccination. Alternatively, White and Medley [13] modelled the increase and decrease of the immunity level and its influence on the pathogen spread in a simple deterministic model, and presented the simulated proportion of hosts according to a continuous immunity level at equilibrium. No model has been developed to study the transient behaviour of the dynamics of the host number according to their immunity level. Such a model would give valuable insight into the interpretation of observed data in terms of, for example, possible time since introduction of infection. This is particularly important for infections that exist in relatively small, discontinuous sub-populations for which local extinction and re-introduction are important processes.
30
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
In this paper, we are interested in modelling the dynamics of immunity levels (i.e. of the number of individuals according to their immunity level) in a cattle population. The deterministic models proposed by White and Medley [13] and by Glass and Grenfell [11] were defined for a large population. However, cattle populations on farms are typically small and tightly controlled. For example, in the United Kingdom, the average cattle herd size is below 150 animals. Consequently, a stochastic approach is required. Decisions to sell or cull are taken by farmers to maintain the size within a small range. The majority of models for pathogen spread do not take into account this controlled size, although the influence on the prevalence of the modelling of the herd size variation in a small range has been previously discussed [14]. The objective of this paper is to represent, for a controlled livestock population, the dynamics of the distribution of individual immunity levels over time by stochastic transmission modelling and to discuss the possible uses of such a model. A general model for diseases with long-lasting immunity is described in Section 2. The structure is a stochastic SIR model with host demography assuming K categories of immunity levels. As 2 or 3 cut-off values are proposed for interpretation of laboratory test results, for example, for ELISA antibody tests, then particular cases are presented in Sections 3 (a model with K = 2 categories of immunity levels: undetectable, detectable) and 4 (a model with K = 3 categories of immunity levels: undetectable, low detectable, high detectable). The numerical examples given in this paper concerned a pathogen inducing a longlasting protection against the disease but with an immunity level which can become undetectable, such as infectious bovine rhinotracheitis (IBR) [15]. Finally, the strengths and weaknesses and potential uses of the model are discussed in Section 5.
2. General model 2.1. Description Let z be the immunity level and t P 0 the time. Because the result of interest is the number of individuals classified by immunity level, and eventually we wish to relate these models to data, the variable z is modelled as a discrete variable. Let K be the number of immunity level categories considered (numbered from 1 to K). Let S(t, z), I(t, z), and R(t, z) be the random variables considered to represent the number of individuals according to their immunity level P categories z at time t in the classes susceptible, infective and immune, respectively, and N ðtÞ ¼ Kz¼1 ðSðt; zÞ þ Iðt; zÞ þ Rðt; zÞÞ the total number of individuals at time t. These variables are non-negative integer valued. Let {X(t) = ({S(t, z)}z=1,. . .,K, {I(t, z)}z=1,. . .,K, {R(t, z)}z=1,. . .,K)}t be the process representing the distribution of the model states considered in this paper. No vertical transmission of the pathogen (from a mother to its offspring during pregnancy) is modelled, nor do we consider pathogen introduction (after the initial condition where at least one infectious animal is introduced at t = 0). The transitions due to each event are represented in Fig. 1 in association with the immunity level, and in Fig. 2 for the particular case where K = 2. The transition parameters for the model are given in Table 1. In our model, susceptible individuals, S, are protected temporarily by passive immunity (maternal antibodies). These individuals are typically born with an immunity level in the non-zero categories (categories numbered higher than 1) which then decreases over time/age. For an
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43 Birth with maternal antibodies
Transmission
31
Re-exposure Recovery
Kth Immunity level categories
2nd 1st Birth without maternal antibodies
S
I
R
Health statuses
Fig. 1. A schematic diagram representing the possible variation of immunity level according to health statuses (S: susceptible, I: infective, R: immune) and the relations with pathogen transmission, recovery and re-exposure. The vertical thick line showed the possible range of immunity levels at time of movement between health statuses.
Immunity level categories
2
S(.,2) kS (2)
1
S(.,1)
I(.,2) k I (1)
λS (.,1, I (.,1), I (.,2), N (.))
I(.,1)
γ (2)
R(.,2) kR (2)
λRe−exp (., z, I (.,1), I (.,2), N (.))
R(.,1)
Fig. 2. Diagrammatic representation of a simple SIR model incorporating immunity levels divided into two categories (undetectable and detectable immunity level). The model’s states are susceptible with undetectable maternal immunity level S(., 1), susceptible with detectable maternal antibody level S(., 2), infective with undetectable immunity level I(., 1), infective with detectable immunity level I(., 2), immune (recovered) with undetectable immunity level R(., 1), and immune (recovered) with detectable immunity level R(., 2). The variable N(.) represents the total size of the population. The transition rates are described Table 1.
individual, the immunity level is assumed to increase when the individual is in the infective status, I, and to decrease when the individual is in the recovered status R. Let kS(z), kI(z) and kR(z) be the parameter-dependent functions representing the rate of decrease of the immunity level for individuals in state S(t, z), the rate of increase of the immunity level for individuals in state I(t, z), and the rate of decrease of the immunity level for individuals in state R(t, z), respectively. Due to boundary conditions, we naturally assume kS(1) = 0, kI(K) = 0 and kR(1) = 0. It is assumed that susceptible individuals with a very high passive immunity level, i.e. in the state S(t, K), are protected against an infection. The transition rate from the state S(t, z) to the state I(t, z), with z < K, should be defined according to the infectivity and the susceptibility of the concerned individuals, both of which vary with immunity. The pathogen transmission is
32
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
Table 1 Possible transitions according to events occurrence in Dt unit of time Event
Transition
Infection of a susceptible (z < K)
S(t, z) ! S(t, z) 1 I(t, z) ! I(t, z) + 1 I(t, z) ! I(t, z) 1 R(t, z) ! R(t, z) + 1 R(t, z) ! R(t, z) 1 R(t, K) ! R(t, K) + 1 S(t, 1) ! S(t, 1) + 1 S(t, z) ! S(t, z) + 1 S(t, z) ! S(t, z) 1 I(t, z) ! I(t, z) 1 R(t, z) ! R(t, z) 1 S(t, z) ! S(t, z) 1 S(t, z 1) ! S(t, z 1) + 1 I(t, z) ! I(t, z) 1 I(t, z + 1) ! I(t,z + 1) + 1 R(t, z) ! R(t, z) 1 R(t, z 1) ! R(t, z 1) + 1
Recovery of an infective (z > 1) Re-exposure of an immune (z < K) Birth of a susceptible without maternal antibodies Birth of a susceptible with maternal antibodies (z > 1) Exit of a susceptible Exit of an infective Exit of an immune Decrease of the immunity level (z > 1) Increase of the immunity level (z < K) Decrease of the immunity level (z > 1)
Transition rate PK bS pðyÞIðt;yÞ y¼1 qS ðzÞSðt; zÞ N ðtÞ c(z)I(t, z) PK bR
y¼1
pðyÞIðt;yÞ
N ðtÞ
qR ðzÞRðt; zÞ
hSS(t, 1) + hII(t, 1) + hRR(t, 1) hII(t, z) + hRR(t, z) lSS(t, z) lII(t, z) lRR(t, z) kS(z)S(t, z) kI(z)I(t, z) kR(z)R(t, z)
assumed to decrease as the immunity level increases. The per-capita rate of infection at time t, kS(t, z, {I(t, z)}z=1,. . .,K, N(t)) is given by P bS Ky¼1 pðyÞIðt; yÞ qS ðzÞ; kS ðt; z; fIðt; zÞgz¼1;...;K ; N ðtÞÞ ¼ N ðtÞ where bS is the infection transmission coefficient, p(z) the relative infectiousness of infective individuals with an immunity level z, z 2 {1, . . . , K}, and qS(z) the relative susceptibility of susceptible individuals with (maternal) immunity level z, z 2 {1, . . . , K}. The functions p(z) and qS(z) decrease over z and we define low immunity values, p(z = 1) = 1 and qS(z = 1) = 1, and high immunity values p(z = K) < 1 and qS(z = K) = 0. After infection, the immunity level of an individual increases during the infective status. It is assumed that animals at time of recovery have an immunity level in the non-zero categories (categories numbered higher than 1). Let c(z) be the recovery rate from the state I(t, z) to the state R(t, z) with z > 1. Due to boundary conditions, we naturally assume c(K) > 0. Once recovered, the immunity level decreases. In case of a re-exposure, the immunity level increases again (boosted immunity). In this case, the immunity level would be assumed to reach the highest possible immunity level category. The transition rate, from the immunity level category z to the Kth category in the immune status R, depends on the infectivity and the susceptibility to re-exposure of the corresponding individuals, both of which vary with immunity. The corresponding per-capita rate, kRe-exp(t, z, {I(t, z)}z=1,. . .,K, N(t)) (re-exposure rate), is given by P bR Ky¼1 pðyÞIðt; yÞ qR ðzÞ; kRe-exp ðt; z; fIðt; zÞgz¼1;...;K ; N ðtÞÞ ¼ N ðtÞ
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
33
where bR is the ‘transmission’ coefficient for re-exposure, p(z) the relative infectiousness of infective individuals with an immunity level z, z 2 {1, . . . , K}, and qR(z) the relative susceptibility of immune animals with an immunity level z, z 2 {1, . . . , K}. Susceptibility to re-exposure decreases over z and we define the lowest and highest immunity levels as qR(z = 1) = 1 and qR(z = K) = 0, respectively. The demography of the population is represented by a simple birth and death process. The birth and the exit rates may depend on the health status and are denoted hS, hI, hR and lS, lI, lR, respectively. Susceptible individuals with maternal antibodies (immunity level categories higher than 1) are assumed to be too young to be pregnant (maternal antibodies are typically lost by 6 months and no pregnancy begun before this age for cattle). Other susceptible individuals are assumed to give birth to susceptible individuals. Non-susceptible individuals are assumed to give birth to susceptible individuals with maternal antibodies which induce the same immunity level as themselves. No specific exit due to the disease is considered (i.e. disease induced neither death nor culling) so all individuals have the same exit rate from the population (lS = lI = lR). 2.2. Analytical formulation Let a = ({sz}z=1,. . .,K, {iz}z=1,. . .,K, {rz}z=1,. . .,K) be one possible value of the continuous time Markov process X(t) = ({S(t, z)}z=1,. . .,K, {I(t, z)}z=1,. . .,K, {R(t, z)}z=1,. . .,K). The state space of {X(t)}tP0 is D = {a = ({sz}z=1,. . .,K, {iz}z=1,. . .,K, {rz}z=1,. . .,K) 2 Z3K, sz, iz, rz P 0 "z = 1, . . . , K}. Let Pa(t) = P(({S(t, z)}z=1,. . .,K, {I(t, z)}z=1,. . .,K, {R(t, z)}z=1,. . .,K) = a) be the probability at time t that the population is in state a 2 D. Assuming no births and no exits, the size of the state space is finite and given by d ¼ C 3K1 N þ3K1 where N is the constant size of the population. To find the exact equilibrium distribution of the model without demography, the eigenvector of the transition rate matrix (d · d matrix) should be obtained, which is difficult in our case. For a very small population of 10 individuals and considering K = 3 categories of immunity levels, we must find an eigenvector of a 43 758 · 43 758 matrix which is a non-trivial task. For the population sizes we are considering, stochastic extinction of the pathogen is likely to be an important process in determining transmission dynamics and control. When the infection is at the quasi-stationary distribution (stationary distribution conditioned on non-extinction), the tranP sition to extinction occurs with a rate equal to lI Qi1 ¼1 þ Kz¼2 ðlI þ cðzÞÞQiz ¼1 where Qiz ¼1 denotes the quasi-stationary probability that there is exactly 1 infective individual in the immunity level z with no other infective individual in the population. As the rate of such transitions is constant over time (conditional upon non-extinction), the distribution of the time to extinction of infection sQ is an exponential distribution with an expected value equal to EðsQ Þ ¼
lI Qi1 ¼1 þ
1 . z¼2 ðlI þ cðzÞÞQiz ¼1
PK
This result is a generalisation of that of Na˚sell [16] for a model with an infectious health status containing many internal stages. The quasi-stationary distribution is more complex to obtain and is not defined here.
34
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
The goal of this work is to describe the dynamics of the distribution of individual immunity levels which depends on the transmission of the pathogen. Let Z(t) = (Z(t, 1), . . . , Z(t, K)) be the distribution of individual immunity levels at time t. The distribution of individual immunity levels {Za(z) = sz + iz + rz}z=1,. . .,K can be obtained from a. Let u = (u1, u2, . . . , uK) be a possible value of Z(t) = (Z(t, 1), . . . , Z(t, K)) and Pu(t) = P(Z(t) = u) be the state probability for u. The corresponding state space is DZ = {u = (u1, u2, . . . , uK) 2 ZK, u1, . . . , uK P 0}. The size of the state space is less than the size of D. Without demography, the state space size is given by d z ¼ C K1 N þK1 where N is the constant size of the population. For a very small population of 10 individuals and considering K = 3 categories of immunity levels, dz = 66. the pathogen spread process, the PP Based on the state probabilities {Pa(t)}a2D forP u(t) could be 0 calculated for all u 2 DZ by the relation P u ðtÞ ¼ a=Z a ¼u P a ðtÞ which implies P u ðtÞ ¼ a=Z a ¼u P 0a ðtÞ, although this is a non-trivial task. Another way is to derive the system of differential equations for these probabilities Pu(t) in the same way as for the Kolmogorov forward equations (not presented here). The pathogen per-capita transmission rate, kS(., z, {I(., z)}z=1,. . .,K, N(.)), and the recovery per-capita rate, c(z), are not directly taken into account in the expression of Pu(t). This is not surprising because the corresponding transitions do not induce a movement between different immunity level categories contrary to the re-exposure transition. Nevertheless, they influence the Pa(t) and thus, indirectly, the Pu(t).
3. Model with K = 2 immunity levels In this section, we consider the special case where the immunity levels can be separated into K = 2 categories: undetectable and detectable immunity levels. The transitions (transmission, recovery, increase and decrease of immunity levels) are represented in Fig. 2. As in the previous section, let a = (s1, s2, i1, i2, r1, r2) be one possible value of X(t) = (S(t, 1), S(t, 2), I(t, 1), I(t, 2), R(t, 1), R(t, 2)) and Pa(t) the state probability for a. 3.1. Epidemiological parameters We do not intend to model a particular infection explicitly, so parameter values used in the simulations shown are illustrative and presented in Table 2. Transmission of the pathogen was assumed to be possible only from infected animals with an undetectable immunity level (level 1) to susceptible animal with an undetectable immunity level. The pathogen was assumed to induce a detectable immunity level (level 2) after 15 days on average (consistent with bovine response to an infection by some pathogens). The infected animals with a detectable immunity level were assumed to recover after 10 days on average. The decrease of the immunity level was assumed to be very slow, on average lasting 2000 days (i.e. more than 5 years). On average, the maternal immunity was assumed to be detectable for 100 days. 3.2. Demographic parameters The birth and exit rates are defined such that the population size remains in a small range over time, and based on values for a dairy cattle herd.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
35
Table 2 Transitions parameters used in the model with 2 immunity levels Event
Parameters
Value
Infection of a susceptible
Coefficient of transmission, bS Transmission level, if immunity level is 1, p(1) Transmission level, if immunity level is 2, p(2) Susceptibility level, if immunity level is 1, qS(1) Susceptibility level, if immunity level is 2, qS(2) Recovery rate, if immunity level is 2, c(2) Coefficient of transmission, bR Susceptibility level, if immunity level is 1, qR(1) Susceptibility level, if immunity level is 2, qR(2) Decrease rate, if immunity level is 2, kS(2) Increase rate, if immunity level is 1, kI(1) Decrease rate, if immunity level is 2, kR(2) Birth rate, h Constant rate, l0 Coefficient of herd size, l1
0.8 1 0 1 0 0.2 0.5 1 0 0.01 0.067 0.0005 0.00061 0.00061 0.00005
Recovery of an infective Re-exposure of an immune
Decrease of the maternal immunity level Increase of the immunity level Decrease of the immunity level Birth Exit
For simplicity, the birth and the exit rates are assumed not to be influenced by the disease so that: l = lS = lI = lR and h = hS = hI = hR. No birth rate is defined for animals with detectable maternal antibodies, i.e. in the model states S(., 2). In order to avoid large variation in herd size, the exit rate was modelled as a linear function of the herd size: l = l0 + l1(N Nref) with l1 constant, and Nref the target population size. The mean exit rate l0 was set to be 0.00061 corresponding to a mean lifespan of 4.5 years, and the birth rate was assumed to have the same value. 3.3. Simulation approach Although the reduced model is simpler than the general one, the size of the state space D is infinite. Without demography, assuming a population size of 50 animals, the size of D is equal to C 550þ5 ¼ 417 451 320. Consequently, the stochastic model is simulated using Matlab (version 6.5) with continuous time (including stochastic simulation of next-event times) as in MacKenzie and Bishop [17]. Over time, the number of individuals in each model state is recorded without tracking of each individual. The initial population consists of 49 susceptible individuals and 1 infective individual, all in the immunity level 1. The model was simulated over 5 years. The transition parameter values used are presented in Table 2. 3.4. Model results 3.4.1. Dynamics of the herd size Based on 500 replications, the minimum, median and maximum herd size stayed within small ranges over the 5 simulated years: [38, 50], [49, 50] and [50, 65], respectively.
36
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
3.4.2. Number of animals in each immunity level The evolution of the number of animals in each immunity level over time for two example replications is presented in Fig. 3. Note that the numbers infected are not explicitly shown in this figure – only the consequent changes in serological immunity. In both examples, the infection fades out (becomes extinct), but there is a considerably difference in the size of the preceding epidemic. Generally, the presence of newly infected animals results in an increase in the number of animals with a high level of immunity, and the lack of sustained, onward transmission results in an increase in low levels of immunity. 3.4.3. Estimation of the time elapsed since the introduction of infection Considering all the replications, the probability over time {Pu(t)}tP0 for each possible distribution of the individual immunity levels u 2 DZ can be obtained. The probability of observing a particular u 2 DZ over time can be used to assess the time since introduction of infection. For example, consider a herd with between 40% and 44% of individuals in the lowest immunity level, i.e. for a herd of 50 animals, this is 20, 21 or 22 individuals with the lowest immunity level. Fig. 4 shows the probability of a herd being in this state over time since introduction. There are two clear modes, so that if this state is observed, it is most likely that the pathogen would have been introduced either a few days ago (A), or between 0.5 and 2 years ago (B).
Fig. 3. Evolution of the number of animals by immunity level over time in two replications: (- - -) undetectable immunity level and (3) detectable immunity level.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
37
Fig. 4. Probability over time to be in the distribution of individual immunity levels corresponding to an observed herd with around 42% of the herd (between 40% and 44%) in the immunity level category 1 (5000 replications).
4. Model with K = 3 immunity levels In this section, we consider the special case where the immunity levels can be separated into K = 3 categories: undetectable, low detectable, and high detectable immunity levels. The transitions (transmission, recovery, increase and decrease of immunity levels) are represented in Fig. 5. As in Section 3, let a = (s1, s2, s3, i1, i2, i3, r1, r2, r3) be one possible value of X(t) = (S(t, 1), S(t, 2), S(t, 3), I(t, 1), I(t, 2), I(t, 3), R(t, 1), R(t, 2), R(t, 3)) and Pa(t) the state probability for a. 4.1. Parameters and simulation As previously, the parameter values used are illustrative. They are presented in Table 3. The demographic parameters are the same as Section 3 (Table 2). Transmission of the pathogen was assumed to be possible from infected animals, with the lowest and intermediate (with a reduced rate) immunity levels, but not from the highest level. Animals with maternal immunity are susceptible with the lowest and intermediate (with a reduced rate) immunity levels but not highest immunity. Other immunological dynamics and demographic parameter values are the same as Section 3.
38
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43 Immunity level categories 3
kI (2)
kS (3)
2
S(.,2)
λS (.,2,{I (., z )}z =1..3 , N (.))
kS (2)
1
S(.,1)
γ (3)
I(.,3)
S(.,3)
I(.,2)
kR (3)
γ (2)
k I (1)
λS (.,1, {I (., z )}z =1..3 , N (.))
R(.,3)
I(.,1)
λRe −exp (., z ,{I (., z )}z=1..3 , N (.))
R(.,2) kR (2)
R(.,1)
Fig. 5. Diagrammatic representation of a simple SIR model incorporating immunity levels divided into 3 categories (undetectable, low detectable and high detectable immunity level). The model’s states are susceptible with undetectable maternal immunity level S(., 1), susceptible with low detectable maternal antibody level S(., 2), susceptible with high detectable maternal antibody level S(., 3), infective with undetectable immunity level I(., 1), infective with low detectable immunity level I(., 2), infective with high detectable immunity level I(., 3), immune (recovered) with undetectable immunity level R(., 1), immune (recovered) with low detectable immunity level R(., 2) and immune (recovered) with high detectable immunity level R(., 3). The variable N(.) represents the total size of the population. The transition rates are described Table 1.
Table 3 Transitions parameters used in the model with 3 immunity levels Event
Parameters
Value
Infection of a susceptible
Coefficient of transmission, bS Transmission level, if immunity level is 1, p(1) Transmission level, if immunity level is 2, p(2) Transmission level, if immunity level is 3, p(3) Susceptibility level, if immunity level is 1, qS(1) Susceptibility level, if immunity level is 2, qS(2) Susceptibility level, if immunity level is 3, qS(3) Recovery rate, if immunity level is 2, c(2) Recovery rate, if immunity level is 3, c(3) Coefficient of transmission, bR Susceptibility level, if immunity level is 1, qR(1) Susceptibility level, if immunity level is 2, qR(2) Susceptibility level, if immunity level is 3, qR(3) Decrease rate, if immunity level is 3, kS(3) Decrease rate, if immunity level is 2, kS(2) Increase rate, if immunity level is 1, kI(1) Increase rate, if immunity level is 2, kI(2) Decrease rate, if immunity level is 3, kR(3) Decrease rate, if immunity level is 2, kR(2)
0.8 1 0.1 0 1 0.5 0 0 0.2 0.5 1 0.5 0 0.02 0.02 0.2 0.067 0.001 0.001
Recovery of an infective Re-exposure of an immune
Decrease of the maternal immunity level Increase of the immunity level Decrease of the immunity level
The stochastic model was simulated in the same way as the model with K = 2. The initial population consists of 49 susceptible animals and 1 infective animal, all in the immunity level category 1. The model was simulated over 5 years.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
39
4.2. Model results 4.2.1. Estimation of the time elapsed since the introduction of infection Fig. 6 shows the probability of observing a particular herd state as a function of time since introduction (as Fig. 4). The particular state considered is between 40% and 44% of animals with undetectable immunity and between 18% and 22% with low detectable immunity. As previously the pathogen may have been introduced either a few days ago (A), or between 1 and 3.5 years ago (B). The inclusion of the additional immune state has the effect of reducing the overall probability of occurrence (compare vertical axis scales), and potentially increasing the precision with which epidemiological quantities can be derived. Note that we have not attempted to standardise the epidemics in Figs. 4 and 6 (e.g. the basic reproduction number is not equal), so that the results are not directly comparable. 4.2.2. Distribution of immunity levels For each immunity level category the probability of occurrence over time can be represented independently (Fig. 7). For the undetectable and the high detectable immunity level categories, the distributions of the number of animals are bimodal during the 3.5 years after the introduction of infection (Fig. 7(a) and (c)). Immediately after introduction, either very few animals or nearly
Fig. 6. Probability over time to be in the distribution of individual immunity levels corresponding to an observed herd with around 42% (between 40% and 44%) and around 20% (between 18% and 22%) of the herd in the immunity level category 1 and 2, respectively (5000 replications).
40
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
Fig. 7. Probability over time of the numbers of animals in the immunity level (a) undetectable (z = 1), (b) low detectable (z = 2), (c) high detectable (z = 3) (500 replications).
all animals are in these extreme categories corresponding to two main dynamics of pathogen spread: early extinction and large epidemic. For the low detectable immunity level category, the bimodal distribution develops more than half a year after the introduction of infection (Fig. 7(b)). This is because the intermediate immune state is derived from the high state, i.e. infected individuals move rapidly to high levels of immunity which then subsequently wanes much more slowly. At the end of the simulation period, the 3 distributions are unimodal, indicating that the distribution of individual immunity levels at this time cannot inform estimates of past infection whatever the total epidemic size that occurred. With these parameters, extinction from the herd within three years is highly likely. Note that the population size changes over time in these simulations, contributing to the variability in the outcome distributions.
5. Discussion The model described in this paper provides insights into the dynamics of the distribution of individual immunity levels at the population level. We believe this to be the first model to consider
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
41
the stochastic dynamics of transmission and immunity. The transmission dynamic model and the basic framework could be developed for specific pathogens by inclusion of, for example, a latency period, vertical transmission, and death due to infection. The subdivision of the health status (S, I and R) according to the immunity level can be seen as the internal staging within each status, so that the residence time in each health status does not follow an exponential distribution. It would also be straight forward to extend the simulation to more immunity categories as may be determined by considering serological data (e.g. undetectable, low, intermediate and high titres). We intend to confront these models with data to assess their usefulness in terms of understanding transient dynamics of infection within herds. Although the basic transmission model is simple, the introduction of the immunity level categories and the transitions between them induce a high complexity, limiting the analytical study of model properties. No exact solution of the stochastic model was obtained, even when the number of immune categories is reduced to 2. In stochastic models with demography but without immunity level categories, only approximations to describe long-term behaviour have been presented (e.g. [16,18]). Such approximations could be sought for the models proposed here, but this is a non-trivial task, and our interests principally concern the transient rather than the long-term behaviour of the pathogen spread and of the distribution of the immunity levels after pathogen introduction. Furthermore, typical livestock (cattle) populations are of such a size that the system of Kolmogorov forward equations will often be unsolvable, so that simulation of the stochastic model is required. Based on a number of replications, the dynamics of the number of individuals in each immunity level can be assessed independently of each other (Fig. 7). A statistical descriptive method to sum up the results of the replications in term of joint distribution of the 3 immunity level categories over time would be useful and should be identified. The model is defined at herd level. The individuals are neither identified nor separated into age groups. Nevertheless, the distribution of the immunity levels according to age groups is particularly informative on the time since introduction of infection since rates of transmission typically vary with age, and age provides a maximum time since infection. Inclusion of age would allow the differentiation between a past pathogen spread and a low pathogen spread: in the first case, some young individuals will have antibodies, whereas in the second case, no young individual will have antibodies. Additionally, farmed livestock are typically divided into groups which, to some extent reflect age, but which reflect production status (e.g. calves, heifers and milking cows and dry cows – i.e. those between lactations). This structuring is used in surveillance, where first screens might include all individuals in the population, and only subsets of young individuals are tested thereafter. To model such data, another approach is needed. Individual-based approach has been used by Glass and Grenfell [11], and Cohen et al. [19]. This approach would increase the complexity of the model. At least, a structure of the population in age groups could be introduced into the model. It would allow informed decisions based on the evolution of the dynamic of the individual immunity levels in the corresponding controlled age groups. In our models, 2 and 3 immunity level categories were considered, corresponding to results of quantitative antibody tests defined with one and two cut-off values, respectively. For defining the possible time periods of pathogen introduction, we assumed that the antibody test used is a gold standard. Nevertheless, available tests often have a specificity and a sensitivity which are less than 100%. To be able to interpret observed data, sensitivity and specificity of the test should be taken into account in the model.
42
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
We have shown that this approach can be used to identify periods of time since pathogen introduction from the set of probabilities of observing the distribution of antibody titre u at time t for all time t {P(ujt)}t>0 (Figs. 4 and 6). The probability P(tju), i.e. the probability of time since introduction given the distribution u can be obtained from: P ðtjuÞ ¼
P ðujtÞP ðtÞ . P ðuÞ
Here, we may take P(t) to be a constant, corresponding to prior ignorance (through note that this is an improper prior) P and the probability, P(u), of observing the distribution of antibody titre u could be obtained as t P ðuktÞP ðtÞ. As the simulation times are finite, P(u) cannot be estimated directly from simulated data, and the use of a Bayesian approach to obtain P(tju) should be further investigated. The precision of simulation results depends on the number of simulations and the uncertainty of parameter values. Sensitivity analyses, a verification and a validation of any model are needed before applying it to specific situations [20]. For verification, the precision of particular outcomes should be related to the number of replications; in particular the results in Figs. 4 and 6 need to be considered in the light of the number of replications used. For validation, available distributions of antibody level for which the date of the pathogen introduction is known could be used to confront the possible introduction periods identified by simulation to the known introduction date. The model can be used to give insight into the observed distribution of the immunity levels, e.g. identifying possible period of times since pathogen introduction. In this paper, only one pathogen introduction was assumed and one population size was considered. If the population size, the routes of the pathogen introduction (e.g. purchase, contact with neighbourhood cattle), and the frequency of pathogen reintroductions (e.g. of purchases of infected animals) are known for an observed distribution, the model can be simulated with the corresponding parameter values. The population size can often be obtained. Nevertheless, identifying the routes of pathogen introductions and the frequencies of pathogen reintroductions is much more difficult. Without such information, a sample of replications from different possible scenarios of simulations with different routes of pathogen introduction and different frequencies of pathogen reintroductions should be run to have a better interpretation of the observed distribution. Acknowledgments The authors thank L.E. Green, L.J. White, and M.J. Keeling for their comments. G.F. Medley acknowledges support from BBSRC (BBS/B/04854). A.-F. Viet was supported by the Institut National de la Recherche Agronomique (France) during the research visit. References [1] H. Houe, Epidemiology of bovine viral diarrhea virus, Vet. Clin. North Am. Food Anim. Pract. 11 (1995) 521. [2] D.J. Nokes, E.A. Okiro, M. Ngama, L.J. White, R. Ochola, P.D. Scott, P.A. Cane, G.F. Medley, RSV epidemiology in a birth cohort in Kilifi District, Kenya: infection in the first year of life, J. Infect. Dis. 190 (2004) 1828.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
43
[3] O. Bahri, M. Ben Halima, M. Ben Ghorbal, K. Dali, Z. Arrouji, S. Khammassi, H. Triki, A. Slim, Measles surveillance and control in Tunisis: 1979–2000, Vaccine 21 (2003) 440. [4] S.H. Maloo, W. Thorpe, G. Kioo, P. Ngumi, G.J. Rowlands, B.D. Perry, Seroprevalences of vector-transmitted infections of small-holder dairy cattle in coastal Kenya, Prev. Vet. Med. 52 (2001) 1. [5] B.T. Grenfell, R.M. Anderson, The estimation of age-related rates of infection from case notifications and serological data, J. Hyg. Camb. 95 (1985) 419. [6] A.E. Ades, D.J. Nokes, Modeling age-specific and time-specific incidence from seroprevalence – toxoplasmosis, Am. J. Epidemiol. 137 (1993) 1022. [7] M.J. Cox, R.S. Azevedo, P.A. Cane, E. Massad, G.F. Medley, Seroepidemiological study of respiratory syncytial virus in Sao Paulo State, Brazil. J. Med. Virol. 55 (1998) 234. [8] N.J. Gay, Analysis of serological surveys using mixture models: application to a survey of parvovirus B19, Stat. Med. 15 (1996) 1567. [9] N.J. Gay, A.J. Vyse, F. Enquselassie, W. Nigatu, D.J. Nokes, Improving sensitivity of oral fluid testing in IgG prevalence studies: application of mixture models to a rubella antibody survey, Epidemiol. Infect. 130 (2003) 285. [10] P.F.M Teunis, O.G. Van Der Heijden, H.E. De Melker, J.F.P. Schellekens, F.G.A. Versteegh, M.E.E. Kretzschmar, Kinetics of the IgG antibody response to pertussis toxin after infection with B. pertussis, Epidemiol. Infect. 129 (2002) 479. [11] K. Glass, B.T. Grenfell, Antibody dynamics in childhood diseases: waning and boosting of immunity and the impact of vaccination, J. Theor. Biol. 221 (2003) 121. [12] M.E.J. Woolhouse, D.T. Haydon, A. Pearson, R.P. Kitching, Failure of vaccination to prevent outbreaks of footand-mouth disease, Epidemiol. Infect. 116 (1996) 363. [13] L.J. White, G.F. Medley, Microparasite population dynamics and continuous immunity, Proc. R. Soc. Lond. B 265 (1998) 1977. [14] D. Clancy, N.P. French, A stochastic model for disease transmission in a managed herd, motivated by Neospora caninum amongst dairy cattle, Math. Biosci. 170 (2001) 113. [15] F. Boelaert, P. Biront, B. Soumare, M. Dispas, E. Vanopdenbosch, J.P. Vermeersch, A. Raskin, J. Dufey, D. Berkvens, P. Kerkhofs, Prevalence of bovine herpesvirus-1 in the Belgian cattle population, Prev. Vet. Med. 45 (2000) 285. [16] I. Na˚sell, Stochastic models of some endemic infections, Math. Biosci. 179 (2002) 1. [17] K. MacKenzie, S.C. Bishop, Developing stochastic models to quantify the dynamics of infectious diseases in domestic livestock, J. Anim. Sci. 79 (2001) 2047. [18] H. Andersson, T. Britton, Stochastic epidemics in dynamic populations: quasi-stationary and extinction, J. Math. Biol. 41 (2000) 559. [19] C. Cohen, M. Artois, D. Pontier, A discrete-event computer model of feline herpes virus within cat populations, Prev. Vet. Med. 45 (2000) 163. [20] J.P.C. Kleijnen, Verification and validation of simulation models, Eur. J. Oper. Res. 82 (1995) 145.
Stochastic dynamics of immunity in small populations: A general framework A.-F. Viet
a,b,*
, G.F. Medley
a
a
b
Ecology and Epidemiology Group, Department of Biological Sciences, University of Warwick, Coventry CV4 7AL, United Kingdom Unit of Animal Health Management, Veterinary School – INRA, BP 40706, 44307 Nantes Cedex 03, France Received 9 August 2004; received in revised form 8 December 2005; accepted 22 December 2005 Available online 8 February 2006
Abstract Assessment of immunological status is a powerful tool in the surveillance and control of infectious pathogens in livestock and human populations. The distribution of immunity levels in the population provides information on time and age dependent transmission. A stochastic model is developed for a livestock population which relates the dynamics of the distribution of immunity levels at the population level to those of pathogen transmission. A general model with K immunity level categories is first proposed, taking into account the increase of the immunity level due to an infection or a re-exposure, the decrease of the immunity level with time since infection or exposure, and the effect of immunity level on the susceptibility and the infectivity of individuals. Numerical results are presented in the particular cases with K = 2 and K = 3 immunity level categories. We demonstrate that for a given distribution of the immunity levels at the population level, the model can be used to identify quantities such as most likely periods of time since introduction of infection. We discuss this approach in relation to analysis of serological data. 2006 Elsevier Inc. All rights reserved. Keywords: Epidemic; Immunity; Stochastic model; Population; Transmission dynamics; Seroepidemiology
*
Corresponding author. Address: Unit of Animal Health Management, Veterinary School – INRA, BP 40706, 44307 Nantes Cedex 03, France. Tel.: +33 (0)240 687 806; fax: +33 (0)240 687 768. E-mail address: [email protected] (A.-F. Viet). 0025-5564/$ - see front matter 2006 Elsevier Inc. All rights reserved. doi:10.1016/j.mbs.2005.12.027
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
29
1. Introduction The rational control of an infectious disease depends on data to assess the incidence and prevalence. Such data usually consist of the health status at time t of the population based either (1) on the record of disease occurrence at individual level before t, or (2) on the presence of an indicator of a previous infection at individual level at time t. The first approach requires that the disease is recognisable and accurately recorded for each individual. However, most infections are asymptomatic or show unspecific signs, for example bovine-viral diarrhoea virus (BVDV) in cattle [1], and respiratory syncytial virus (RSV) in humans [2]. Consequently, specific diagnosis frequently relies on serology (i.e. finding specific antibodies to the infection within an individual). The second approach relies on surveying the population for pathogen-specific antibody levels as an indicator of a previous infection. Such data give a measure of the extent of previous infection and, indirectly, of the prevalence at the time of the survey. For examples of the use of such serology see Bahri et al. [3] and Maloo et al. [4]. Cross-sectional and longitudinal surveys of pathogen-specific antibody levels provide valuable data for assessment of infection dynamics, especially when combined with information on age to provide age serological profiles (e.g. [5–7]). The immunity level of an individual is determined by measurement of the concentration of specific antibodies to the pathogen in a sample of serum or, more recently, oral fluid (‘saliva’). For most clinical needs, it is sufficient to categorise individuals into two classes: positive (and protected) and negative (and susceptible to infection). As assessment of antibody concentration is quantitative (titre), this requires that the original continuous measurement is reclassified into a binary score based on a predefined or relative cut-off value. However, it has recently been realised that the distribution of titres itself contains valuable information [8,9]. Because titre generally increases during or immediately after an infection or a re-exposure and decreases with time since exposure, the antibody level potentially contains information on the time since infection of the individual. Thus, distribution of specific individual antibody levels contains indirect information on the transmission dynamics of the pathogen. For example, if the pathogen has been present, but is now absent, a few individuals would have a high immunity level and most individuals a low level, whereas, with the pathogen present, most individuals would have a high immunity level. Previous work has considered the relation between pathogen transmission and the distribution of the immunity level at the population level. Some researchers aimed at interpreting such observed distributions in terms of pathogen spread or of time since infection of the population (e.g. [1,10]). In contrast, Glass and Grenfell [11] and Woolhouse et al. [12] introduced antibody levels to consider waning protection in childhood disease and against infection with foot and mouth disease virus, respectively, and the impact of vaccination. Alternatively, White and Medley [13] modelled the increase and decrease of the immunity level and its influence on the pathogen spread in a simple deterministic model, and presented the simulated proportion of hosts according to a continuous immunity level at equilibrium. No model has been developed to study the transient behaviour of the dynamics of the host number according to their immunity level. Such a model would give valuable insight into the interpretation of observed data in terms of, for example, possible time since introduction of infection. This is particularly important for infections that exist in relatively small, discontinuous sub-populations for which local extinction and re-introduction are important processes.
30
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
In this paper, we are interested in modelling the dynamics of immunity levels (i.e. of the number of individuals according to their immunity level) in a cattle population. The deterministic models proposed by White and Medley [13] and by Glass and Grenfell [11] were defined for a large population. However, cattle populations on farms are typically small and tightly controlled. For example, in the United Kingdom, the average cattle herd size is below 150 animals. Consequently, a stochastic approach is required. Decisions to sell or cull are taken by farmers to maintain the size within a small range. The majority of models for pathogen spread do not take into account this controlled size, although the influence on the prevalence of the modelling of the herd size variation in a small range has been previously discussed [14]. The objective of this paper is to represent, for a controlled livestock population, the dynamics of the distribution of individual immunity levels over time by stochastic transmission modelling and to discuss the possible uses of such a model. A general model for diseases with long-lasting immunity is described in Section 2. The structure is a stochastic SIR model with host demography assuming K categories of immunity levels. As 2 or 3 cut-off values are proposed for interpretation of laboratory test results, for example, for ELISA antibody tests, then particular cases are presented in Sections 3 (a model with K = 2 categories of immunity levels: undetectable, detectable) and 4 (a model with K = 3 categories of immunity levels: undetectable, low detectable, high detectable). The numerical examples given in this paper concerned a pathogen inducing a longlasting protection against the disease but with an immunity level which can become undetectable, such as infectious bovine rhinotracheitis (IBR) [15]. Finally, the strengths and weaknesses and potential uses of the model are discussed in Section 5.
2. General model 2.1. Description Let z be the immunity level and t P 0 the time. Because the result of interest is the number of individuals classified by immunity level, and eventually we wish to relate these models to data, the variable z is modelled as a discrete variable. Let K be the number of immunity level categories considered (numbered from 1 to K). Let S(t, z), I(t, z), and R(t, z) be the random variables considered to represent the number of individuals according to their immunity level P categories z at time t in the classes susceptible, infective and immune, respectively, and N ðtÞ ¼ Kz¼1 ðSðt; zÞ þ Iðt; zÞ þ Rðt; zÞÞ the total number of individuals at time t. These variables are non-negative integer valued. Let {X(t) = ({S(t, z)}z=1,. . .,K, {I(t, z)}z=1,. . .,K, {R(t, z)}z=1,. . .,K)}t be the process representing the distribution of the model states considered in this paper. No vertical transmission of the pathogen (from a mother to its offspring during pregnancy) is modelled, nor do we consider pathogen introduction (after the initial condition where at least one infectious animal is introduced at t = 0). The transitions due to each event are represented in Fig. 1 in association with the immunity level, and in Fig. 2 for the particular case where K = 2. The transition parameters for the model are given in Table 1. In our model, susceptible individuals, S, are protected temporarily by passive immunity (maternal antibodies). These individuals are typically born with an immunity level in the non-zero categories (categories numbered higher than 1) which then decreases over time/age. For an
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43 Birth with maternal antibodies
Transmission
31
Re-exposure Recovery
Kth Immunity level categories
2nd 1st Birth without maternal antibodies
S
I
R
Health statuses
Fig. 1. A schematic diagram representing the possible variation of immunity level according to health statuses (S: susceptible, I: infective, R: immune) and the relations with pathogen transmission, recovery and re-exposure. The vertical thick line showed the possible range of immunity levels at time of movement between health statuses.
Immunity level categories
2
S(.,2) kS (2)
1
S(.,1)
I(.,2) k I (1)
λS (.,1, I (.,1), I (.,2), N (.))
I(.,1)
γ (2)
R(.,2) kR (2)
λRe−exp (., z, I (.,1), I (.,2), N (.))
R(.,1)
Fig. 2. Diagrammatic representation of a simple SIR model incorporating immunity levels divided into two categories (undetectable and detectable immunity level). The model’s states are susceptible with undetectable maternal immunity level S(., 1), susceptible with detectable maternal antibody level S(., 2), infective with undetectable immunity level I(., 1), infective with detectable immunity level I(., 2), immune (recovered) with undetectable immunity level R(., 1), and immune (recovered) with detectable immunity level R(., 2). The variable N(.) represents the total size of the population. The transition rates are described Table 1.
individual, the immunity level is assumed to increase when the individual is in the infective status, I, and to decrease when the individual is in the recovered status R. Let kS(z), kI(z) and kR(z) be the parameter-dependent functions representing the rate of decrease of the immunity level for individuals in state S(t, z), the rate of increase of the immunity level for individuals in state I(t, z), and the rate of decrease of the immunity level for individuals in state R(t, z), respectively. Due to boundary conditions, we naturally assume kS(1) = 0, kI(K) = 0 and kR(1) = 0. It is assumed that susceptible individuals with a very high passive immunity level, i.e. in the state S(t, K), are protected against an infection. The transition rate from the state S(t, z) to the state I(t, z), with z < K, should be defined according to the infectivity and the susceptibility of the concerned individuals, both of which vary with immunity. The pathogen transmission is
32
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
Table 1 Possible transitions according to events occurrence in Dt unit of time Event
Transition
Infection of a susceptible (z < K)
S(t, z) ! S(t, z) 1 I(t, z) ! I(t, z) + 1 I(t, z) ! I(t, z) 1 R(t, z) ! R(t, z) + 1 R(t, z) ! R(t, z) 1 R(t, K) ! R(t, K) + 1 S(t, 1) ! S(t, 1) + 1 S(t, z) ! S(t, z) + 1 S(t, z) ! S(t, z) 1 I(t, z) ! I(t, z) 1 R(t, z) ! R(t, z) 1 S(t, z) ! S(t, z) 1 S(t, z 1) ! S(t, z 1) + 1 I(t, z) ! I(t, z) 1 I(t, z + 1) ! I(t,z + 1) + 1 R(t, z) ! R(t, z) 1 R(t, z 1) ! R(t, z 1) + 1
Recovery of an infective (z > 1) Re-exposure of an immune (z < K) Birth of a susceptible without maternal antibodies Birth of a susceptible with maternal antibodies (z > 1) Exit of a susceptible Exit of an infective Exit of an immune Decrease of the immunity level (z > 1) Increase of the immunity level (z < K) Decrease of the immunity level (z > 1)
Transition rate PK bS pðyÞIðt;yÞ y¼1 qS ðzÞSðt; zÞ N ðtÞ c(z)I(t, z) PK bR
y¼1
pðyÞIðt;yÞ
N ðtÞ
qR ðzÞRðt; zÞ
hSS(t, 1) + hII(t, 1) + hRR(t, 1) hII(t, z) + hRR(t, z) lSS(t, z) lII(t, z) lRR(t, z) kS(z)S(t, z) kI(z)I(t, z) kR(z)R(t, z)
assumed to decrease as the immunity level increases. The per-capita rate of infection at time t, kS(t, z, {I(t, z)}z=1,. . .,K, N(t)) is given by P bS Ky¼1 pðyÞIðt; yÞ qS ðzÞ; kS ðt; z; fIðt; zÞgz¼1;...;K ; N ðtÞÞ ¼ N ðtÞ where bS is the infection transmission coefficient, p(z) the relative infectiousness of infective individuals with an immunity level z, z 2 {1, . . . , K}, and qS(z) the relative susceptibility of susceptible individuals with (maternal) immunity level z, z 2 {1, . . . , K}. The functions p(z) and qS(z) decrease over z and we define low immunity values, p(z = 1) = 1 and qS(z = 1) = 1, and high immunity values p(z = K) < 1 and qS(z = K) = 0. After infection, the immunity level of an individual increases during the infective status. It is assumed that animals at time of recovery have an immunity level in the non-zero categories (categories numbered higher than 1). Let c(z) be the recovery rate from the state I(t, z) to the state R(t, z) with z > 1. Due to boundary conditions, we naturally assume c(K) > 0. Once recovered, the immunity level decreases. In case of a re-exposure, the immunity level increases again (boosted immunity). In this case, the immunity level would be assumed to reach the highest possible immunity level category. The transition rate, from the immunity level category z to the Kth category in the immune status R, depends on the infectivity and the susceptibility to re-exposure of the corresponding individuals, both of which vary with immunity. The corresponding per-capita rate, kRe-exp(t, z, {I(t, z)}z=1,. . .,K, N(t)) (re-exposure rate), is given by P bR Ky¼1 pðyÞIðt; yÞ qR ðzÞ; kRe-exp ðt; z; fIðt; zÞgz¼1;...;K ; N ðtÞÞ ¼ N ðtÞ
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
33
where bR is the ‘transmission’ coefficient for re-exposure, p(z) the relative infectiousness of infective individuals with an immunity level z, z 2 {1, . . . , K}, and qR(z) the relative susceptibility of immune animals with an immunity level z, z 2 {1, . . . , K}. Susceptibility to re-exposure decreases over z and we define the lowest and highest immunity levels as qR(z = 1) = 1 and qR(z = K) = 0, respectively. The demography of the population is represented by a simple birth and death process. The birth and the exit rates may depend on the health status and are denoted hS, hI, hR and lS, lI, lR, respectively. Susceptible individuals with maternal antibodies (immunity level categories higher than 1) are assumed to be too young to be pregnant (maternal antibodies are typically lost by 6 months and no pregnancy begun before this age for cattle). Other susceptible individuals are assumed to give birth to susceptible individuals. Non-susceptible individuals are assumed to give birth to susceptible individuals with maternal antibodies which induce the same immunity level as themselves. No specific exit due to the disease is considered (i.e. disease induced neither death nor culling) so all individuals have the same exit rate from the population (lS = lI = lR). 2.2. Analytical formulation Let a = ({sz}z=1,. . .,K, {iz}z=1,. . .,K, {rz}z=1,. . .,K) be one possible value of the continuous time Markov process X(t) = ({S(t, z)}z=1,. . .,K, {I(t, z)}z=1,. . .,K, {R(t, z)}z=1,. . .,K). The state space of {X(t)}tP0 is D = {a = ({sz}z=1,. . .,K, {iz}z=1,. . .,K, {rz}z=1,. . .,K) 2 Z3K, sz, iz, rz P 0 "z = 1, . . . , K}. Let Pa(t) = P(({S(t, z)}z=1,. . .,K, {I(t, z)}z=1,. . .,K, {R(t, z)}z=1,. . .,K) = a) be the probability at time t that the population is in state a 2 D. Assuming no births and no exits, the size of the state space is finite and given by d ¼ C 3K1 N þ3K1 where N is the constant size of the population. To find the exact equilibrium distribution of the model without demography, the eigenvector of the transition rate matrix (d · d matrix) should be obtained, which is difficult in our case. For a very small population of 10 individuals and considering K = 3 categories of immunity levels, we must find an eigenvector of a 43 758 · 43 758 matrix which is a non-trivial task. For the population sizes we are considering, stochastic extinction of the pathogen is likely to be an important process in determining transmission dynamics and control. When the infection is at the quasi-stationary distribution (stationary distribution conditioned on non-extinction), the tranP sition to extinction occurs with a rate equal to lI Qi1 ¼1 þ Kz¼2 ðlI þ cðzÞÞQiz ¼1 where Qiz ¼1 denotes the quasi-stationary probability that there is exactly 1 infective individual in the immunity level z with no other infective individual in the population. As the rate of such transitions is constant over time (conditional upon non-extinction), the distribution of the time to extinction of infection sQ is an exponential distribution with an expected value equal to EðsQ Þ ¼
lI Qi1 ¼1 þ
1 . z¼2 ðlI þ cðzÞÞQiz ¼1
PK
This result is a generalisation of that of Na˚sell [16] for a model with an infectious health status containing many internal stages. The quasi-stationary distribution is more complex to obtain and is not defined here.
34
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
The goal of this work is to describe the dynamics of the distribution of individual immunity levels which depends on the transmission of the pathogen. Let Z(t) = (Z(t, 1), . . . , Z(t, K)) be the distribution of individual immunity levels at time t. The distribution of individual immunity levels {Za(z) = sz + iz + rz}z=1,. . .,K can be obtained from a. Let u = (u1, u2, . . . , uK) be a possible value of Z(t) = (Z(t, 1), . . . , Z(t, K)) and Pu(t) = P(Z(t) = u) be the state probability for u. The corresponding state space is DZ = {u = (u1, u2, . . . , uK) 2 ZK, u1, . . . , uK P 0}. The size of the state space is less than the size of D. Without demography, the state space size is given by d z ¼ C K1 N þK1 where N is the constant size of the population. For a very small population of 10 individuals and considering K = 3 categories of immunity levels, dz = 66. the pathogen spread process, the PP Based on the state probabilities {Pa(t)}a2D forP u(t) could be 0 calculated for all u 2 DZ by the relation P u ðtÞ ¼ a=Z a ¼u P a ðtÞ which implies P u ðtÞ ¼ a=Z a ¼u P 0a ðtÞ, although this is a non-trivial task. Another way is to derive the system of differential equations for these probabilities Pu(t) in the same way as for the Kolmogorov forward equations (not presented here). The pathogen per-capita transmission rate, kS(., z, {I(., z)}z=1,. . .,K, N(.)), and the recovery per-capita rate, c(z), are not directly taken into account in the expression of Pu(t). This is not surprising because the corresponding transitions do not induce a movement between different immunity level categories contrary to the re-exposure transition. Nevertheless, they influence the Pa(t) and thus, indirectly, the Pu(t).
3. Model with K = 2 immunity levels In this section, we consider the special case where the immunity levels can be separated into K = 2 categories: undetectable and detectable immunity levels. The transitions (transmission, recovery, increase and decrease of immunity levels) are represented in Fig. 2. As in the previous section, let a = (s1, s2, i1, i2, r1, r2) be one possible value of X(t) = (S(t, 1), S(t, 2), I(t, 1), I(t, 2), R(t, 1), R(t, 2)) and Pa(t) the state probability for a. 3.1. Epidemiological parameters We do not intend to model a particular infection explicitly, so parameter values used in the simulations shown are illustrative and presented in Table 2. Transmission of the pathogen was assumed to be possible only from infected animals with an undetectable immunity level (level 1) to susceptible animal with an undetectable immunity level. The pathogen was assumed to induce a detectable immunity level (level 2) after 15 days on average (consistent with bovine response to an infection by some pathogens). The infected animals with a detectable immunity level were assumed to recover after 10 days on average. The decrease of the immunity level was assumed to be very slow, on average lasting 2000 days (i.e. more than 5 years). On average, the maternal immunity was assumed to be detectable for 100 days. 3.2. Demographic parameters The birth and exit rates are defined such that the population size remains in a small range over time, and based on values for a dairy cattle herd.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
35
Table 2 Transitions parameters used in the model with 2 immunity levels Event
Parameters
Value
Infection of a susceptible
Coefficient of transmission, bS Transmission level, if immunity level is 1, p(1) Transmission level, if immunity level is 2, p(2) Susceptibility level, if immunity level is 1, qS(1) Susceptibility level, if immunity level is 2, qS(2) Recovery rate, if immunity level is 2, c(2) Coefficient of transmission, bR Susceptibility level, if immunity level is 1, qR(1) Susceptibility level, if immunity level is 2, qR(2) Decrease rate, if immunity level is 2, kS(2) Increase rate, if immunity level is 1, kI(1) Decrease rate, if immunity level is 2, kR(2) Birth rate, h Constant rate, l0 Coefficient of herd size, l1
0.8 1 0 1 0 0.2 0.5 1 0 0.01 0.067 0.0005 0.00061 0.00061 0.00005
Recovery of an infective Re-exposure of an immune
Decrease of the maternal immunity level Increase of the immunity level Decrease of the immunity level Birth Exit
For simplicity, the birth and the exit rates are assumed not to be influenced by the disease so that: l = lS = lI = lR and h = hS = hI = hR. No birth rate is defined for animals with detectable maternal antibodies, i.e. in the model states S(., 2). In order to avoid large variation in herd size, the exit rate was modelled as a linear function of the herd size: l = l0 + l1(N Nref) with l1 constant, and Nref the target population size. The mean exit rate l0 was set to be 0.00061 corresponding to a mean lifespan of 4.5 years, and the birth rate was assumed to have the same value. 3.3. Simulation approach Although the reduced model is simpler than the general one, the size of the state space D is infinite. Without demography, assuming a population size of 50 animals, the size of D is equal to C 550þ5 ¼ 417 451 320. Consequently, the stochastic model is simulated using Matlab (version 6.5) with continuous time (including stochastic simulation of next-event times) as in MacKenzie and Bishop [17]. Over time, the number of individuals in each model state is recorded without tracking of each individual. The initial population consists of 49 susceptible individuals and 1 infective individual, all in the immunity level 1. The model was simulated over 5 years. The transition parameter values used are presented in Table 2. 3.4. Model results 3.4.1. Dynamics of the herd size Based on 500 replications, the minimum, median and maximum herd size stayed within small ranges over the 5 simulated years: [38, 50], [49, 50] and [50, 65], respectively.
36
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
3.4.2. Number of animals in each immunity level The evolution of the number of animals in each immunity level over time for two example replications is presented in Fig. 3. Note that the numbers infected are not explicitly shown in this figure – only the consequent changes in serological immunity. In both examples, the infection fades out (becomes extinct), but there is a considerably difference in the size of the preceding epidemic. Generally, the presence of newly infected animals results in an increase in the number of animals with a high level of immunity, and the lack of sustained, onward transmission results in an increase in low levels of immunity. 3.4.3. Estimation of the time elapsed since the introduction of infection Considering all the replications, the probability over time {Pu(t)}tP0 for each possible distribution of the individual immunity levels u 2 DZ can be obtained. The probability of observing a particular u 2 DZ over time can be used to assess the time since introduction of infection. For example, consider a herd with between 40% and 44% of individuals in the lowest immunity level, i.e. for a herd of 50 animals, this is 20, 21 or 22 individuals with the lowest immunity level. Fig. 4 shows the probability of a herd being in this state over time since introduction. There are two clear modes, so that if this state is observed, it is most likely that the pathogen would have been introduced either a few days ago (A), or between 0.5 and 2 years ago (B).
Fig. 3. Evolution of the number of animals by immunity level over time in two replications: (- - -) undetectable immunity level and (3) detectable immunity level.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
37
Fig. 4. Probability over time to be in the distribution of individual immunity levels corresponding to an observed herd with around 42% of the herd (between 40% and 44%) in the immunity level category 1 (5000 replications).
4. Model with K = 3 immunity levels In this section, we consider the special case where the immunity levels can be separated into K = 3 categories: undetectable, low detectable, and high detectable immunity levels. The transitions (transmission, recovery, increase and decrease of immunity levels) are represented in Fig. 5. As in Section 3, let a = (s1, s2, s3, i1, i2, i3, r1, r2, r3) be one possible value of X(t) = (S(t, 1), S(t, 2), S(t, 3), I(t, 1), I(t, 2), I(t, 3), R(t, 1), R(t, 2), R(t, 3)) and Pa(t) the state probability for a. 4.1. Parameters and simulation As previously, the parameter values used are illustrative. They are presented in Table 3. The demographic parameters are the same as Section 3 (Table 2). Transmission of the pathogen was assumed to be possible from infected animals, with the lowest and intermediate (with a reduced rate) immunity levels, but not from the highest level. Animals with maternal immunity are susceptible with the lowest and intermediate (with a reduced rate) immunity levels but not highest immunity. Other immunological dynamics and demographic parameter values are the same as Section 3.
38
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43 Immunity level categories 3
kI (2)
kS (3)
2
S(.,2)
λS (.,2,{I (., z )}z =1..3 , N (.))
kS (2)
1
S(.,1)
γ (3)
I(.,3)
S(.,3)
I(.,2)
kR (3)
γ (2)
k I (1)
λS (.,1, {I (., z )}z =1..3 , N (.))
R(.,3)
I(.,1)
λRe −exp (., z ,{I (., z )}z=1..3 , N (.))
R(.,2) kR (2)
R(.,1)
Fig. 5. Diagrammatic representation of a simple SIR model incorporating immunity levels divided into 3 categories (undetectable, low detectable and high detectable immunity level). The model’s states are susceptible with undetectable maternal immunity level S(., 1), susceptible with low detectable maternal antibody level S(., 2), susceptible with high detectable maternal antibody level S(., 3), infective with undetectable immunity level I(., 1), infective with low detectable immunity level I(., 2), infective with high detectable immunity level I(., 3), immune (recovered) with undetectable immunity level R(., 1), immune (recovered) with low detectable immunity level R(., 2) and immune (recovered) with high detectable immunity level R(., 3). The variable N(.) represents the total size of the population. The transition rates are described Table 1.
Table 3 Transitions parameters used in the model with 3 immunity levels Event
Parameters
Value
Infection of a susceptible
Coefficient of transmission, bS Transmission level, if immunity level is 1, p(1) Transmission level, if immunity level is 2, p(2) Transmission level, if immunity level is 3, p(3) Susceptibility level, if immunity level is 1, qS(1) Susceptibility level, if immunity level is 2, qS(2) Susceptibility level, if immunity level is 3, qS(3) Recovery rate, if immunity level is 2, c(2) Recovery rate, if immunity level is 3, c(3) Coefficient of transmission, bR Susceptibility level, if immunity level is 1, qR(1) Susceptibility level, if immunity level is 2, qR(2) Susceptibility level, if immunity level is 3, qR(3) Decrease rate, if immunity level is 3, kS(3) Decrease rate, if immunity level is 2, kS(2) Increase rate, if immunity level is 1, kI(1) Increase rate, if immunity level is 2, kI(2) Decrease rate, if immunity level is 3, kR(3) Decrease rate, if immunity level is 2, kR(2)
0.8 1 0.1 0 1 0.5 0 0 0.2 0.5 1 0.5 0 0.02 0.02 0.2 0.067 0.001 0.001
Recovery of an infective Re-exposure of an immune
Decrease of the maternal immunity level Increase of the immunity level Decrease of the immunity level
The stochastic model was simulated in the same way as the model with K = 2. The initial population consists of 49 susceptible animals and 1 infective animal, all in the immunity level category 1. The model was simulated over 5 years.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
39
4.2. Model results 4.2.1. Estimation of the time elapsed since the introduction of infection Fig. 6 shows the probability of observing a particular herd state as a function of time since introduction (as Fig. 4). The particular state considered is between 40% and 44% of animals with undetectable immunity and between 18% and 22% with low detectable immunity. As previously the pathogen may have been introduced either a few days ago (A), or between 1 and 3.5 years ago (B). The inclusion of the additional immune state has the effect of reducing the overall probability of occurrence (compare vertical axis scales), and potentially increasing the precision with which epidemiological quantities can be derived. Note that we have not attempted to standardise the epidemics in Figs. 4 and 6 (e.g. the basic reproduction number is not equal), so that the results are not directly comparable. 4.2.2. Distribution of immunity levels For each immunity level category the probability of occurrence over time can be represented independently (Fig. 7). For the undetectable and the high detectable immunity level categories, the distributions of the number of animals are bimodal during the 3.5 years after the introduction of infection (Fig. 7(a) and (c)). Immediately after introduction, either very few animals or nearly
Fig. 6. Probability over time to be in the distribution of individual immunity levels corresponding to an observed herd with around 42% (between 40% and 44%) and around 20% (between 18% and 22%) of the herd in the immunity level category 1 and 2, respectively (5000 replications).
40
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
Fig. 7. Probability over time of the numbers of animals in the immunity level (a) undetectable (z = 1), (b) low detectable (z = 2), (c) high detectable (z = 3) (500 replications).
all animals are in these extreme categories corresponding to two main dynamics of pathogen spread: early extinction and large epidemic. For the low detectable immunity level category, the bimodal distribution develops more than half a year after the introduction of infection (Fig. 7(b)). This is because the intermediate immune state is derived from the high state, i.e. infected individuals move rapidly to high levels of immunity which then subsequently wanes much more slowly. At the end of the simulation period, the 3 distributions are unimodal, indicating that the distribution of individual immunity levels at this time cannot inform estimates of past infection whatever the total epidemic size that occurred. With these parameters, extinction from the herd within three years is highly likely. Note that the population size changes over time in these simulations, contributing to the variability in the outcome distributions.
5. Discussion The model described in this paper provides insights into the dynamics of the distribution of individual immunity levels at the population level. We believe this to be the first model to consider
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
41
the stochastic dynamics of transmission and immunity. The transmission dynamic model and the basic framework could be developed for specific pathogens by inclusion of, for example, a latency period, vertical transmission, and death due to infection. The subdivision of the health status (S, I and R) according to the immunity level can be seen as the internal staging within each status, so that the residence time in each health status does not follow an exponential distribution. It would also be straight forward to extend the simulation to more immunity categories as may be determined by considering serological data (e.g. undetectable, low, intermediate and high titres). We intend to confront these models with data to assess their usefulness in terms of understanding transient dynamics of infection within herds. Although the basic transmission model is simple, the introduction of the immunity level categories and the transitions between them induce a high complexity, limiting the analytical study of model properties. No exact solution of the stochastic model was obtained, even when the number of immune categories is reduced to 2. In stochastic models with demography but without immunity level categories, only approximations to describe long-term behaviour have been presented (e.g. [16,18]). Such approximations could be sought for the models proposed here, but this is a non-trivial task, and our interests principally concern the transient rather than the long-term behaviour of the pathogen spread and of the distribution of the immunity levels after pathogen introduction. Furthermore, typical livestock (cattle) populations are of such a size that the system of Kolmogorov forward equations will often be unsolvable, so that simulation of the stochastic model is required. Based on a number of replications, the dynamics of the number of individuals in each immunity level can be assessed independently of each other (Fig. 7). A statistical descriptive method to sum up the results of the replications in term of joint distribution of the 3 immunity level categories over time would be useful and should be identified. The model is defined at herd level. The individuals are neither identified nor separated into age groups. Nevertheless, the distribution of the immunity levels according to age groups is particularly informative on the time since introduction of infection since rates of transmission typically vary with age, and age provides a maximum time since infection. Inclusion of age would allow the differentiation between a past pathogen spread and a low pathogen spread: in the first case, some young individuals will have antibodies, whereas in the second case, no young individual will have antibodies. Additionally, farmed livestock are typically divided into groups which, to some extent reflect age, but which reflect production status (e.g. calves, heifers and milking cows and dry cows – i.e. those between lactations). This structuring is used in surveillance, where first screens might include all individuals in the population, and only subsets of young individuals are tested thereafter. To model such data, another approach is needed. Individual-based approach has been used by Glass and Grenfell [11], and Cohen et al. [19]. This approach would increase the complexity of the model. At least, a structure of the population in age groups could be introduced into the model. It would allow informed decisions based on the evolution of the dynamic of the individual immunity levels in the corresponding controlled age groups. In our models, 2 and 3 immunity level categories were considered, corresponding to results of quantitative antibody tests defined with one and two cut-off values, respectively. For defining the possible time periods of pathogen introduction, we assumed that the antibody test used is a gold standard. Nevertheless, available tests often have a specificity and a sensitivity which are less than 100%. To be able to interpret observed data, sensitivity and specificity of the test should be taken into account in the model.
42
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
We have shown that this approach can be used to identify periods of time since pathogen introduction from the set of probabilities of observing the distribution of antibody titre u at time t for all time t {P(ujt)}t>0 (Figs. 4 and 6). The probability P(tju), i.e. the probability of time since introduction given the distribution u can be obtained from: P ðtjuÞ ¼
P ðujtÞP ðtÞ . P ðuÞ
Here, we may take P(t) to be a constant, corresponding to prior ignorance (through note that this is an improper prior) P and the probability, P(u), of observing the distribution of antibody titre u could be obtained as t P ðuktÞP ðtÞ. As the simulation times are finite, P(u) cannot be estimated directly from simulated data, and the use of a Bayesian approach to obtain P(tju) should be further investigated. The precision of simulation results depends on the number of simulations and the uncertainty of parameter values. Sensitivity analyses, a verification and a validation of any model are needed before applying it to specific situations [20]. For verification, the precision of particular outcomes should be related to the number of replications; in particular the results in Figs. 4 and 6 need to be considered in the light of the number of replications used. For validation, available distributions of antibody level for which the date of the pathogen introduction is known could be used to confront the possible introduction periods identified by simulation to the known introduction date. The model can be used to give insight into the observed distribution of the immunity levels, e.g. identifying possible period of times since pathogen introduction. In this paper, only one pathogen introduction was assumed and one population size was considered. If the population size, the routes of the pathogen introduction (e.g. purchase, contact with neighbourhood cattle), and the frequency of pathogen reintroductions (e.g. of purchases of infected animals) are known for an observed distribution, the model can be simulated with the corresponding parameter values. The population size can often be obtained. Nevertheless, identifying the routes of pathogen introductions and the frequencies of pathogen reintroductions is much more difficult. Without such information, a sample of replications from different possible scenarios of simulations with different routes of pathogen introduction and different frequencies of pathogen reintroductions should be run to have a better interpretation of the observed distribution. Acknowledgments The authors thank L.E. Green, L.J. White, and M.J. Keeling for their comments. G.F. Medley acknowledges support from BBSRC (BBS/B/04854). A.-F. Viet was supported by the Institut National de la Recherche Agronomique (France) during the research visit. References [1] H. Houe, Epidemiology of bovine viral diarrhea virus, Vet. Clin. North Am. Food Anim. Pract. 11 (1995) 521. [2] D.J. Nokes, E.A. Okiro, M. Ngama, L.J. White, R. Ochola, P.D. Scott, P.A. Cane, G.F. Medley, RSV epidemiology in a birth cohort in Kilifi District, Kenya: infection in the first year of life, J. Infect. Dis. 190 (2004) 1828.
A.-F. Viet, G.F. Medley / Mathematical Biosciences 200 (2006) 28–43
43
[3] O. Bahri, M. Ben Halima, M. Ben Ghorbal, K. Dali, Z. Arrouji, S. Khammassi, H. Triki, A. Slim, Measles surveillance and control in Tunisis: 1979–2000, Vaccine 21 (2003) 440. [4] S.H. Maloo, W. Thorpe, G. Kioo, P. Ngumi, G.J. Rowlands, B.D. Perry, Seroprevalences of vector-transmitted infections of small-holder dairy cattle in coastal Kenya, Prev. Vet. Med. 52 (2001) 1. [5] B.T. Grenfell, R.M. Anderson, The estimation of age-related rates of infection from case notifications and serological data, J. Hyg. Camb. 95 (1985) 419. [6] A.E. Ades, D.J. Nokes, Modeling age-specific and time-specific incidence from seroprevalence – toxoplasmosis, Am. J. Epidemiol. 137 (1993) 1022. [7] M.J. Cox, R.S. Azevedo, P.A. Cane, E. Massad, G.F. Medley, Seroepidemiological study of respiratory syncytial virus in Sao Paulo State, Brazil. J. Med. Virol. 55 (1998) 234. [8] N.J. Gay, Analysis of serological surveys using mixture models: application to a survey of parvovirus B19, Stat. Med. 15 (1996) 1567. [9] N.J. Gay, A.J. Vyse, F. Enquselassie, W. Nigatu, D.J. Nokes, Improving sensitivity of oral fluid testing in IgG prevalence studies: application of mixture models to a rubella antibody survey, Epidemiol. Infect. 130 (2003) 285. [10] P.F.M Teunis, O.G. Van Der Heijden, H.E. De Melker, J.F.P. Schellekens, F.G.A. Versteegh, M.E.E. Kretzschmar, Kinetics of the IgG antibody response to pertussis toxin after infection with B. pertussis, Epidemiol. Infect. 129 (2002) 479. [11] K. Glass, B.T. Grenfell, Antibody dynamics in childhood diseases: waning and boosting of immunity and the impact of vaccination, J. Theor. Biol. 221 (2003) 121. [12] M.E.J. Woolhouse, D.T. Haydon, A. Pearson, R.P. Kitching, Failure of vaccination to prevent outbreaks of footand-mouth disease, Epidemiol. Infect. 116 (1996) 363. [13] L.J. White, G.F. Medley, Microparasite population dynamics and continuous immunity, Proc. R. Soc. Lond. B 265 (1998) 1977. [14] D. Clancy, N.P. French, A stochastic model for disease transmission in a managed herd, motivated by Neospora caninum amongst dairy cattle, Math. Biosci. 170 (2001) 113. [15] F. Boelaert, P. Biront, B. Soumare, M. Dispas, E. Vanopdenbosch, J.P. Vermeersch, A. Raskin, J. Dufey, D. Berkvens, P. Kerkhofs, Prevalence of bovine herpesvirus-1 in the Belgian cattle population, Prev. Vet. Med. 45 (2000) 285. [16] I. Na˚sell, Stochastic models of some endemic infections, Math. Biosci. 179 (2002) 1. [17] K. MacKenzie, S.C. Bishop, Developing stochastic models to quantify the dynamics of infectious diseases in domestic livestock, J. Anim. Sci. 79 (2001) 2047. [18] H. Andersson, T. Britton, Stochastic epidemics in dynamic populations: quasi-stationary and extinction, J. Math. Biol. 41 (2000) 559. [19] C. Cohen, M. Artois, D. Pontier, A discrete-event computer model of feline herpes virus within cat populations, Prev. Vet. Med. 45 (2000) 163. [20] J.P.C. Kleijnen, Verification and validation of simulation models, Eur. J. Oper. Res. 82 (1995) 145.