key: cord-0701836-9wxa6jcl
authors: Tuckwell, Henry C.; Toubiana, Laurent
title: Dynamical modeling of viral spread in spatially distributed populations: stochastic origins of oscillations and density dependence
date: 2006-12-20
journal: Biosystems
DOI: 10.1016/j.biosystems.2006.12.006
sha: 90a5ee465cb9a3b549a37a8e71d98ac5b064acef
doc_id: 701836
cord_uid: 9wxa6jcl

In order to understand the spatio-temporal structure of epidemics beyond that permitted with classical SIR (susceptible-infective-recovered)-type models, a new mathematical model for the spread of a viral disease in a population of spatially distributed hosts is described. The positions of the hosts are randomly generated in a rectangular habitat. Encounters between any pair of individuals are according to a Poisson process with a mean rate that declines exponentially as the distance between them increases. The contact rate allows the mean rates to be set at a certain number of encounters per day on average. The relevant state variables for each individual at any time are given by the solution of a pair of coupled differential equations for the viral load and the quantity of general immune system effectors which reduce the viral load. The parameters describing within-host viral-immune system dynamics are generated randomly to reflect variability across a population. Transmission is assumed to depend on the viral loads in donors and occurs with a probability [Formula: see text]. The initial conditions are such that one randomly chosen individual carries a randomly chosen amount of the virus, whereas the rest of the population is uninfected. Simulations reveal local or whole-population responses. Whole-population disease spread may be in the form of isolated or multiple occurrences, the latter often being approximately periodic. The mechanisms of this oscillatory behaviour are analyzed in terms of several parameters and the distribution of critical points in the host dynamical systems. Increased contact rate, increased probability of transmission and decreased threshold for viral transmission, decreased immune strength and increased viral growth rate all increase the probability of multiple outbreaks and the distribution of the critical points also plays a role.

Viruses constitute some of the most dangerous threats to human health. They may invade a community and spread rapidly amongst its members, sometimes causing a large number of fatalities, possibly on a recurrent basis. In the case of the influenza virus, in 1918-1919 a new strain ("Spanish", type H1N1) swept the globe, causing an estimated 50-100 million deaths (de la Barrera and Reyes-Teran, 2005; Osterholm, 2005; Strauss et al., 2006) . At the present time there are several countries, particularly in Africa, with up to 35% of their populations between the ages of 15 and 50 years infected by human immunodeficiency virus (HIV) . Throughout the world, there are currently approximately 40 million persons infected (Regoes et al., 2005) and already over 16 million deaths have been caused by this virus.

Renewed interest in the dynamical processes involved in the spread of viruses has arisen recently due to the threat of bioterrorist attacks, especially with such viruses as smallpox (Henderson et al., 1999; Kaplan et al., 2002) which was eradicated many years ago. More recently, certain parts of the world were thrown into panic and economic chaos by outbreaks of SARS (severe acute respiratory syndrome) in which certain new highly contagious strains of a type of corona virus were implicated. Global spread was facilitated by airline travel, and the number of cases doubled about every week in certain Asian countries in the early stages of this epidemic (Lipsitch et al., 2003; Anon., 2006b) . At the time of writing (early 2006) Western Europe is, as is the rest of the world, in the grip of the fear that "bird flu" (avian influenza H5N1) will move further west and its virus will mutate to a form which is easily transmitted from human to human, thus raising the threat of a global pandemic in which an estimated number of deaths as high as 180-360 million may occur (de la Barrera and Reyes-Teran, 2005; Osterholm, 2005) with economic costs in the US alone estimated at between 70 and 170 billion dollars.

The classically employed models in the analysis of epidemic data are of the SIR type such as those of Kermack-McKendrick (Bailey, 1975; Kermack and McKendrick, 1927) , in which the time variable is continuous and the model is formulated as a system of ordinary differential equations, or discrete models such as Reed-Frost (Abbey, 1952) , extended recently in Tuckwell and Williams (2007) . Examples of such systems with as many as 10 components applied to analysis of vaccination against carrier borne diseases are summarized in Anderson and May (1991) . The SIR approach, which is based on an assumption of homogeneous mixing of susceptible and infectious individuals, is used to make predictions in relation to both practical and theoretical aspects of epidemics and the spread of disease (see Diekmann and Heesterbeek (2000) , Hethcote (2000) , Pollard (1973) for reviews). However, such gross compartmental models do not allow for the inclusion of either the spatial distributions of the individual members of the population nor the detailed dynamics of the infecting particles (viruses or bacteria).

There have been several attempts to improve on simplified models in recent years (see, for example, the review in Isham (2005) ); for example, populations have been divided into groups or patches such that in each subgroup an SIS or SEIR (susceptibleexposed-infective-recovered) model applies, sometimes with migration between patches (Aparicio et al., 2000; Arrigoni and Pugliese, 2002; Lloyd and Jansen, 2004) and threshold conditions determined (Arino et al., 2005; Arino and Van den Driessche, 2003; Wang and Zhao, 2004) , although stochastic effects were often ignored. Longini and coworkers Patel et al., 2005; Weycker et al., 2005) have included stochastic mechanisms and age groups in models for the spread of influenza in communities of total population 10,000, with a view to examining the effects of various vaccination schedules. Furthermore, large scale simulations for populations of many millions (EpiSims) based on realistic social networks (Eubank et al., 2004; Stroud et al., 2006) have highlighted the deficiencies contained in the homogeneous mixing assumption of the classical models.

In addition to assigning spatial locations to the hosts, we will incorporate statistical distributions for the parameters describing the invading pathogen's dynamics, thus allowing for differences in the immune properties of various individuals. Such factors must have a considerable influence on the evolution of a disease. In this context, the idea of "superspreaders" has been discussed in relationship to SARS (Lipsitch et al., 2003) , though the emphasis was on contact rates rather than host viral-immune dynamical parameters. In previous publications (Tuckwell et al., , 2001 we have adopted this new approach to the spread of a disease using a set of dynamical variables which describe the disease states of individuals (animals or humans) through the size of the actual viral populations they harbor and components of the immune system. It is recognized that often there are many locations of viral populations in a single host, but these are considered together. We have considered a network of n > 1 hosts within each of which a viral population may establish itself and such that virus particles may be transmitted from any diseased host to another host. The model we will consider places emphasis on viruses as the disease-causing agent, but our approach may be also applied to bacterial diseases.

The immune system response to viral or other infection is too complicated to attempt to model exactly in every detail. Generally, there are agents in the immune system which, when stimulated by antigens, give rise to host defence mechanisms which are of several types, but whose net action is to eliminate the virus. Active defense mechanisms are lumped into effectors (Anderson and May, 1991, Ch. 3; Bocharov and Romanyukha, 1994) so that a simplified two-dimensional system of effectors and virus may be considered in each individual. Since the immune responses are specific to any particular virus, the use of a system of two differential equations to model the evolution of the virus-effector populations is a simplification but this approach is expected to contain the essential dynamics of the biological processes without rendering the system of equations unwieldy with a very large parameter set. Our approach can easily incorporate the "standard model" of viral dynamics in which there are three components, including uninfected T-cells, infected T-cells and virus particles (Neumann et al., 1998) , or other similar models (Wodarz et al., 2002) . However, this would require even more parameters so here we concentrate on the simplest comprehensive model and reserve extensions to later work. An article using our approach but with the inclusion of immunological memory has recently appeared (Kostova, 2007) .

Thus, when there is no interaction between hosts, a host virus population v(t) and a host effector population a(t) at time t evolve in a simplified picture according to the coupled equations (Anderson and May, 1991; Tuckwell et al., 1998; Tuckwell and Wan, 2000) dv

where r is the intrinsic viral per capita growth rate, γ the clearance rate of virus per effector per virus, λ the rate of production of effectors (for example T-cells), μ the per capita removal (death) rate of effectors, and determines the rate of production of effectors, in response to infection, per unit amount of virus. In general, when interaction between hosts is allowed via transmission of virus particles from one to the other, we may write equations such as the following for the the viral load in the ith host, v i (t) at time t, and the effector population, a i (t), in the ith host

where f i and g i describe the dynamics of the viral and effector populations in the ith individual. Here the matrix element β ij , which may depend on several factors, describes the strength of transmission from host j to host i, and the function F, which in the simplest cases is linear, describes the overall effect of transmission to a host from the rest of the population. One usually takes β ii = 0 for all i. In France, (and many other countries) data on influenza, and other diseases, are obtained from a preassigned group of general practitioners (500 of the approximately 50,000 in the whole country) who record the number of new cases of patients with "flu-like symptoms" and as such can be used to form a rough estimate of the actual numbers of cases. Records of the number of recorded cases of influenza (Anon., 2006a) show major roughly periodic outbreaks, which are probably attributable in part to seasonal effects (Lloyd and Jansen, 2004) , although the precise mechanisms involved are still not understood. One of the questions we wish to address in the present work concerns the nature and causes of periodic outbreaks. We wish to enquire whether "epidemics" may arise due to spatio-temporal stochastic processes in the transmission process. For example, if per chance a spatial cluster of diseased hosts occurs, there may ensue a local outbreak of considerable magnitude which may propagate to the rest of the population. We have already considered the effects of noise in individual responses and noted a broad-based stochastic resonance as the noise parameter changes in nonlinear models with threshold . We wish to investigate here and in subsequent studies the factors, both demographic and dynamical, which influence the frequency, magnitude and duration of various outbreaks.

The model employed here differs from those employed previously in several aspects. A stochastic representation is employed for the locations of members of the host population over a two-dimensional habitat, A. Members have "home" locations which are interpreted as their average positions, from which they may make excursions and encounter other members of the community. The coordinates (X i , Y i ), i = 1, 2, . . . , n, of the "home" locations of the n individuals comprizing the population are assumed to be independent and random over A. (We adopt the notation that random variables are in capital letters, whereas the values they assume are in lower case. However, not all capitalized variables are random.) In the following analysis, distances may be taken to be in kilometers and we assume that A is a rectangle with sides of length a and b. However, the actual physical dimensions are secondary to the rates of contact between individuals. For the first studies, X i and Y i are taken to be uniformly distributed on (0, a) and (0, b), respectively. See Fig. 1 for a typical random spatial distribution with n = 100 individuals and a = b = 1.

Encounters between individuals are assumed to occur at the event times of n(n−1) 2 Poisson processes N ij , i, j = 1, . . . , n; j > i. The Poisson process N ij = {N ij (t), t ≥ 0} governing encounters between individuals i and j is assumed to have rate parameter given by

is the distance between the ith and jth home locations. The parameter Λ is chosen to be proportional to the number of meetings per unit time. No intra-day variations in interaction rate are incorporated. Poisson processes are the natural choice for random events such as chance meetings (Tuckwell, 1995) . The stochastic differential equations governing the evolution of the viral and effector populations in the ith individual are then given by

Here T ji describes the strength of transmission from the jth to the ith individual if a meeting occurs between them. The term dN ji dt indicates that when an event occurs in the Poisson process N ji (signifying a meeting has occurred between individuals j and i) then, since N ji jumps by unity, the derivative of V i has a delta-function at the time of the meeting so that the viral load V i jumps by the amount of virus transmitted from individual j to individual i.

For simplicity in our first investigations we assume that the transmitted amount of virus from individual j to i is given by

independent of the recipient, where β is a constant and G is a given function. This could also be made a function of the dynamical states of donor and recipient but in order to limit the number of parameters we have chosen β to be same for all pairs of individuals. The form of G chosen in the simulations reported below is that of a step function

so that an infected host is capable of transmitting virus particles only if his viral load is greater than the critical level v crit . However, in the simulations, transmission only occurs with probability p trans . This means that if an infected host encounters another host, which occurs when an event happens in the Poisson process governing their meetings, then if the infected host's level of virus is above the threshold v crit for transmission, then virus will be transmitted with probability p trans . The transmission rate can also be endowed with an upper limit such as provided by a logistic function of V j , but here in the first instance for simplicity we assume (4A) is valid. For our first investigations of this mathematical model, as explained in Section 1, the functions f and g describing the intra-host dynamics of the effectorviral populations, are chosen as in our previous articles (Anderson and May, 1991; Tuckwell et al., 1998) so that in the absence of interactions:

In simulations we often modify (5) slightly by making the growth term logistic to prevent the occurrence of unrealistically large viral loads-see below.

The system (5) and (6) has equilibria at

or generically at P 1 and P 2 . In order to take into account variability across the host population, the parameters r i , γ i , λ i , μ i and i are randomly distributed. All these parameters are nonnegative but in accordance with the most natural choice are given normal distributions (central limit theorem, Tuckwell, 1995) with means and standard deviations such that negative values have practically zero chance to occur-and in fact never did occur in any of the trials we have performed. The standard values of the means and standard deviations of these five parameters are shown in Table 1 , their units being as follows: λ i , effectors days −1 ; μ i , days −1 ; r i , days −1 ; γ i , effectors −1 days −1 ; and i , virions −1 days −1 . These values are based in part on (Anderson and May, 1991 ) and on our previous work (Tuckwell et al., , 2001 and are such that when they assume their mean values, the solutions of (5) and (6) in response to an initial small increase in viral load, are such that the viral load rises to a maximum after a few days and then declines to negligible values after about 10 days, while the effector level rises and remains at a higher level (immune memory); that is there are no oscillations. The actual magnitudes of the maximum viral loads are known in some cases where this is the form of the viral population growth curve, e.g. influenza in man (Bocharov and Romanyukha, 1994) and lymphocytic choriomeningitis virus in mice (Bocharov, 1998) . In these references very complicated models (10 to over 40 differential equations) for viral and immune system dynamics have been given but we have employed a reduced system which captures the essence of the very large systems. The key variables (viral load and effector numbers) can be easily rescaled if required to actual biological values (Bocharov, 1998; Bocharov and Romanyukha, 1994; Neumann et al., 1998; Sidorenko, 2005) . For example, viral densities in influenza are of order 10 8 virus particles per millilitre of infected tissue (Sidorenko, 2005) and in lymphocytic choriomeningitis virus in mice, between 10 7 and 10 9 per millilitre (Bocharov, 1998) . In Fig.  2 is shown an example of the randomly distributed Table 1 . In the top part is shown the distribution for the P 1 's for each individual in the population and in the bottom part that for the P 2 's. Note that some of the latter are at unphysical values as they have negative values of the viral load -see Tuckwell and Wan (2000) for a detailed discussion -but this has no effect as the dynamical variables are constrained to have nonnegative values. Here the population is n = 100. positions of the critical points P 1 and P 2 in the (a,v)plane.

In Fig. 3 are shown histograms of the numbers of critical points of various types associated with each individual in the population, which in this example has a size of 100. There are two possibilities for P 1 : unstable saddle point or stable node whereas for P 2 there are Fig. 3 . For the critical points whose coordinates were given in the previous figure we here show the numbers of each kind of critical point. Upper histogram for P 1 and lower histogram for P 2 . three possibilities: unstable saddle, stable node and stable spiral point. Note that if λ i γ i > μ i r i the point P 2,i occurs at negative v and hence is not at a biologically relevant value. There is then just one meaningful critical point P 1 on the a-axis and this is an asymptotically stable node-see Tuckwell and Wan (2000) for a detailed discussion.

In Table 2 are summarized the parameters which complement those of Table 1 . Some of the viral dynamical parameters appear here that are not randomly distributed as are those in Table 1 .

It is remarked that many persons infected with, for example, influenza virus show no symptoms but may still transmit the disease (Regoes et al., 2005) . This is incorporated in our model as v crit < v s . Some remaining demographic parameters are given in Table 3 .

The coordinates (X i , Y i ) of the n individuals in the area of interest are determined using a random number generator. This enables the distance d ij between individual i and j to be determined. An initial distribution of infected individuals is chosen across the population, usually consisting of just one infected individual as depicted in Fig. 1 . This individual is chosen randomly from the population and so has a randomly chosen set of viral dynamical parameters as in Table 1 . The value of the viral load in the initially infected individual is chosen from a uniform distribution on (0, v init ), where, based on , the value of v init was usually (5)) in the absence of transmissions of virus from other infected individuals. All individuals in the population initially have zero effectors: A i (0) = 0 for all i. The rates λ ij , where i = 1, . . . , n, and j = 1, . . . , n, of meetings for each pair of individuals can be determined using Eq. (1). Suppose the time interval on which the spread and evolution of the virus is to be observed is (0, T ] with a given initial condition and with a given time step t. At each time step, a meetings matrix M of order n × n is then generated where M ij (which depends on the time step) is such that

The value of M ij is set at unity if a uniform on (0, 1) random variable is less than λ ij t, in accordance with the Poisson law. Note that M is a symmetric matrix (M ij = M ji ) because in this model, if individual i meets individual j, then individual j meets individual i. Furthermore, care is taken to ensure λ ij t < 1. At each time step the ith individual's effector population A i is first updated according to

if individual i has never been infected,

It is then ascertained if individual i, for i = 1, . . . , n, had any encounters at all with other individuals. If there were no meetings, so that j M ij = 0, then V i evolves according to

Here the parameter v max (constant, non-random) is employed in the growth equation for host viral populations in the absence of viral transmission such that growth is logistic (Bocharov, 1998; Tuckwell et al., 1998) according to the generic equatioṅ

If there is at least one meeting, then the viral load V j (t) in each encountered individual j, j = 1, . . . , n, j = i, is examined to see if it exceeded the critical value v crit . If it did then transmission occurs to individual i with probability p trans with V i being increased by the amount β if transmission was successful. In summary the contribution from individual j to i if they meet is

where U(0, 1) is a random variable uniformly distributed on the interval (0, 1), and T ii (t) = 0 for all i. Thus if individual i encounters anyone he receives the total random viral dosē

However, it is possible thatT i (t) is zero because of the step functions in Eq. (11) in which case V i (t) changes according to Eq. (9).

Broadly speaking, simulations of the above model resulted in one of the following outcomes: response here means that individuals do become infected by any amount of the virus:

(a) a small local response which dies away; (b) a response which involves a large fraction, and possibly all, of the population where individuals become sick (v > v s ) once only, with viral levels returning thereafter to almost zero; (c) as in (b) except that viral levels remain elevated in certain individuals at lower than symptomatic (v < v s ) or infectious (v < v crit ) levels; (d) a repetitive outbreak of disease in which a large fraction, and possibly all, of the population becomes ill on an approximately periodic basis.

An example of a case (d)-type response where there are almost-periodic outbreaks is shown in Fig. 4 , where the mean viral load (in scaled units) across the population is plotted against time and is seen to have approximately periodic peaks. For these results we have used the standard parameters with n = 100 and p trans = 0.2. The spatial distribution of hosts and the distributions of the critical points are those shown in Figs. 1-3 . Recall that only one member of the population is initially infected and that the individual responses are not oscillatory when the viral dynamical parameters take their standard mean values. However, due to the distribution of parameters several individuals may have recurrent outbreaks of diminishing magnitude. However, it was found that recurrent outbreaks in individuals were by themselves not intrinsically sufficient to cause another major outbreak because when the meetings between individuals were stopped at 50 days, there was just the initial outbreak even though all parameters and the spatial distribution were unchanged. In Fig. 5 is shown the temporal evolution of the viral loads in two of the individual hosts for the first outbreak in Fig.  4 . The blue curve, starting at v > 0 at t = 0, depicts the viral load for the initially infected individual, the (other) red curve being for an individual who became infected about 4-5 days later. With a smaller value of p trans = 0.1 and all other parameters the same, with no cut-off of meetings after 50 days, there was only one outbreak.

In order to consider the influence of various parameters, we may classify them into three main types:

(i) those which determine contact rates between hosts, being n, a, b, Λ and α; (ii) those which govern the processes of viral transmission, which are p trans , β and v crit ; and (iii) those which determine the within-host viralimmune system dynamics, being the means and variances of the parameters λ i , μ i , r i , γ i and i , i = 1, . . . , n.

We do not have space to report the dependence of the course of disease on all these in detail, so we discuss variations in a few representative ones. Our simulations were done with software which does not permit the efficient study of very large populations. However we expect that similar principles will hold in the investigations reported here as in larger populations. In later articles we will carry out further investigations of the parameter dependence and incorporate larger population sizes. In all the following results there is assumed to be just one individual at t = 0 who carries a small viral load where both the individual and the viral load are assigned randomly (see Section 3).

In order to obtain an idea of the effects of population density on the multiplicity of outbreaks, 10 trials were first performed with p trans = 0.5 and n = 30, n = 40 and n = 50. Note that since the standard area is one unit, n gives the population density. In these runs, the values of the host dynamical parameters (Table 1) were the same (after being randomly generated) and the spatial distributions of the individuals were also the same (after being randomly generated). It was found that there was some variability in the response type, but often the result was stereotypical. For example, with n = 30, in 8 trials of 10, there were 6 outbreaks in the period t < 250 and in 2 trials there were 5 outbreaks in the same time period. If the viral dynamical parameters in Table 1 are fixed and the spatial distribution is fixed, the only remaining sources of variability are in the meeting times and the success or failure of transmission of the virus in an encounter with an infected individual. For the 10 trials with n = 50, however, there were in every trial almost regular periodic outbreaks.

Including randomly distributed dynamical parameters and random spatial configurations should give a clearer picture of the stochastic nature of the responses. Hence 20 such trials were performed with n = 20, 25, 30 and 35 for 0 ≤ t ≤ 400 days. The number of trials was here relatively small because with the software employed, such runs were very time consuming. In the bar diagrams of Fig. 6 we show distributions, from these simulations, of the multiplicities of the outbreaks for these various values of n. It can be seen that as n increases, the general trend is for the distribution to move to the right where the multiplicity increases. When n = 20 only 10greater than 2, whereas when n = 35 the percentage is 70 and when n = 50 it is 100%.

The hypothesis of an increase in the tendency to multiple outbreaks with increased population density is reinforced by the results given in Fig. 7. Here is plotted the number of times (in 20 trials) that only a single outbreak occurred versus population density from n = 20 to n = 60, with 20 trials for each value of n. These results include the four first column results of Fig. 6 . When n = 20 there were 15 occurrences of just one outbreak, for n = 25 just 11, for n = 30 there were 8, for n = 35 there were 3 and for n = 60, only one in 20 trials gave rise to a single outbreak.

The fact that increasing n leads to a greater chance of multiple outbreaks can be traced to the effect of n on the rates of meetings for any individual which for fixed values of all other parameters are very close to proportional to n. This assertion was verified by performing several simulations with only variations in the parameter Λ-see Section 4.3 for related results. In order to understand which other factors might be primarily responsible for the occurrence of a multiplicity of outbreaks, we examined the nature of the critical points which are of the kind given in the scheme of Fig. 3 . If one or more members of the population have, for the critical point P 2 , asymptotically stable spiral points as opposed to nodes, one might expect that there is a chance that the resulting individual oscillatory activity could manifest itself in multiple outbreaks across the entire population. That is, recurrent epidemics might arise just because one or a few individuals have immune system parameters which lead to their having recurrent illness, which enables them to transmit virus particles to to others. However, it was found that recurrent outbreaks may occur in the population when there are just a few spiral P 2 's or if there are many. 

We now consider the effects of varying the following two parameters:

(a) the probability p trans of a successful transmission of virus from a diseased individual to a host he or she encounters, and, (b) v crit , the value of the viral load required to make a host infectious; that is, capable of transmitting virus to another individual he encounters.

In order to obtain more insight into the stochastic mechanisms involved in disease spread, 20 runs were performed with n = 30 and several values of p trans for each of two sets of host viral dynamical parameters (chosen randomly from the distributions with means and variances given in Table 1 ) and different spatial distributions. We will refer to these two lots of parameter sets and spatial configurations as sets 1 and 2. All other parameters took their standard values in Tables 2 and 3 . For the various values of p trans , the proportion of trials was recorded in which multiple periodic outbreaks, of the kind shown in Fig. 4 , occurred. The proportion of trials thus obtained gives an estimate of the probability of multiple periodic outbreaks. The results are given in Table 4 and indicate the importance of the parameter p trans in influencing whether there will just be a single, possibly major, outbreak, followed by almost no further disease processes, or a sequence of several nearly periodic outbreaks.

For the first set of results, when p trans = 0.2 or 0.25, the chance of recurrent periodic outbreaks is negligible. However, when p trans is increased to 0.3, 50% of trials resulted in recurrent outbreaks, presumably because a new outbreak is easier to trigger with an increased probability of transmission of the virus from infected to other individuals. With p trans = 0.4, 90% of trials gave recurrent outbreaks and when p trans was 0.5, all of the trials performed resulted in sustained oscillatory responses, of the kind depicted in Fig. 4 . For the second set of results periodic outbreaks were harder to elicit, as p trans had to be greater than 0.4 for them to occur at all. Examination of the spatial distributions for the two sets of results showed that for the first set, the initially infected individual was near the center of mass of the population whereas for the second set the initially infected individual was isolated at the periphery. Since the greater the distance between two individuals the smaller their rate of meeting, the position of the initially infected individual may have influenced the probability of multiple periodic outbreaks. (See Sidorenko (2005) for a discussion of related edge effects on the size and duration of epidemics.) It seems from these few results that stochastic mechanisms play a role in determining whether an epidemic consists of one major outbreak or whether there will be oscillations of the viral load across the population.

To illustrate the different kinds of results, Fig. 8 shows the time-courses for the number of sick individuals, defined here as those having a viral load greater than v s (see Table 2 ) in the population for several values of the transmission probability, p trans , which are marked on the figure. Here the parameters and all initial conditions, including the spatial distribution are identical to those employed for set 1 results in Table 4 .

We also investigated the effects of changing the parameter v crit -see Section 2. With n = 30, 20 trials were performed, with randomly chosen dynamical parameters and random spatial configurations. The aver- age numbers of outbreaks with v crit values of 2, 4, 6 and 8 were 2.6, 1.7, 1.1 and 0.9, respectively. This suggests that higher values of v crit tend to make multiple outbreaks less likely. In fact, for a given set of parameters, including the whole spatial distribution, there is a critical value of v crit above which multiple outbreaks tend not to occur.

The parameters which govern the within-host dynamics are the means and variances of the λ i , μ i , r i , γ i and i , i = 1, . . . , n. To ascertain the effects of greater host immunity and of more aggressive virus, we describe results for the effects of varying the means of two of these five parameters; firstly, λ i , which describes the strength of effector production and second, r i , which describes the growth rates of the virus within the hosts. (Analysis given in Tuckwell and Wan (2000) as well as simulation shows that for several parameters their ratio mainly determines the behaviour of solutions, so increases in one has a similar effect to decreases in the another.)

Simulations were performed to see the effects of changes in the mean of λ i , the host intrinsic effector production rates. At the same time contact rates Λ were varied to see the interplay between these parameters. The standard value of the mean E[λ i ] = 0.5 and the values 0.4 and 0.5 were also employed. As the mean E[λ i ] varies, the nature of the distribution of critical points changes, with larger E[λ i ] producing less spiral points. The other parameters in these simulations were standard, with n = 20 and p trans = 0.5. The results based on 20 trials are given in Table 5 .

It can be seen that when the value of E[λ i ] is small, which indicates poor immune status, multiple outbreaks (affecting the whole population) are relatively easy to elicit. In all cases the trend is clear, that larger val- ues of Λ promote multiple outbreaks. Even when the immune status is strong with E[λ i ] = 0.6, multiple outbreaks occurred if the contact rate was high enough. Note that the occurrence or not of multiple outbreaks was not strongly influenced by the number of spiral points for P 2 , indicating a stochastic mechanism.

Simulations were performed with p trans = 0.5 for a population of n = 50 in which the mean of the viral growth rates r i varied from its standard value of 1 to half this value and double this value. All other parameters had their standard values and several runs were performed for 0 ≤ t ≤ 600 days. The results for each trial were similar, probably because the distributions of critical points were similar. Typical results are shown in Fig. 9 . When the mean E[r i ] = 0.5, there is only one initial outbreak, affecting the whole population (top graph). In the case E[r i ] = 1, (the standard value) there are several outbreaks with an interval of a few 100 days between them (middle graph). However, when the mean of r i is 2, as in the bottom plot, the frequency of outbreaks increases greatly. The parameter E[r i ] clearly has a strong influence on the temporal evolution of the disease spread. Mutation to a more virulent form of virus could thus lead to a greater multiplicity of outbreaks.

The model considered in this article captures some of the features of natural populations in that individuals (human or animals), who are potential hosts for disease producing agents such as bacteria or virus particles, are assigned spatial locations and meet other individuals, who may be infected, at random times. It has been assumed that each individual has a home base and his or her interactions with other individuals occur randomly at the event times in Poisson processes whose rates depend on the distances between the home bases. The manner in which locations are chosen could make the model appropriate for certain plant populations. Small-world network effects (Weycker et al., 2005) which may occur in human populations are captured to some degree because any individual may encounter any other individual with a non-zero probability no matter how great the distance between them. As a general rule however, in human and especially animal populations, encounters are on average more likely between individuals the smaller the distance between them, as described by the present model.

The density of the population can be varied and also the contact rates between the individuals. Disease particles are transmitted stochastically from infected to other individuals. The novel feature of our model is that it includes dynamical systems for within-host disease particles (viruses or bacteria) and the host immune response. This has been done in a simple way by using a system of two differential equations which captures the elements of effector-viral or effector-bacterial dynamics. With the use of such models one may distinguish more accurately the spread of various diseases by using appropriate parameters. For example, models for HIV within-host dyamics are well documented (Neumann et al., 1998) and could be used in place of the simple model ( (5) and (6)) and could also be made stochastic (Le Corfec and Tuckwell, 1998; Merrill, 2005) .

As elaborated on in Section 4, the parameters may be classified into three main types: first, those which determine contact rate; second, those which govern the processes of transmission; and those which determine the within-host disease dynamics being the means and variances of the λ i , μ i , r i , γ i and i , i = 1, . . . , n. The number of parameters is thus very large. There are 10 fixed parameters given in Tables 2 and 3, n(n−1) 2 rates of contact, and 10n viral-immune system dynamical parameters-the latter being generated from the 5 means and standard deviations of each of the parameters in Table 1 for each member of the population. Thus each host has an assigned set of parameters so that different hosts will respond to disease agents in a different way. With such a large number of diverse parameters, it is not so simple to find the classically defined basic reproductive ratio R 0 (Anderson and May, 1991; Lloyd and Jansen, 2004) although a more detailed investigation could yield bounds on an estimate as in the large urban simulations (Eubank et al., 2004; Stroud et al., 2006) . We focussed on the occurrence or not of multiple approximately periodic outbreaks and investigated one or more of each class of parameter which seemed to play a role in determining whether single or multiple outbreaks would occur. For the simulations we performed, values of p trans below 0.3 never gave multiple outbreaks (see Table 4 ) whereas values equal to or greater than 0.5 could give multiple outbreaks. The parameter v crit was also investigated and found to exert considerable control over the multiplicity of outbreaks with a welldefined threshold effect. Another property investigated was population density, and since the area of habitat was fixed, it was sufficient to vary the population size n. The main influence of an increase in population density is, if the parameters Λ and α are left unchanged, to increase the contact rates between hosts. It was found that the probability of a single outbreak, as opposed to multiple outbreaks, dropped rapidly as the population density increased. Similarly, if the parameter Λ increased there was an increase in the likelihood of a multiple outbreak. The natures of the population critical points was found to have some influence, but contact rates seemed to be the more important of these factors: if there were a few or many P 2 spiral points multiple outbreaks occurred if Λ was large enough.

Two of the within-host dynamical parameters were varied and found to have an important influence. First the mean of the intrinsic effector production rate, E[λ i ] was varied to just above and below its standard value. Larger E[λ i ] made it harder, but not impossible if Λ increased, to elicit multiple outbreaks, reflecting stronger immune properties of the host population. Smaller values increased the likelihood of recurring illness across the population. Changes in the mean of r i , the withinhost viral production rate had the opposite effects to the corresponding changes in the mean of λ i . In general, the the mean values of the within-host dynamical systems parameters have been based on a restricted set of values, so it remains to be seen what epidemic behaviour other sets of values will induce. In general such parameters, or their analogues in more detailed models, will be disease-specific. This could accommodate the fact that in some diseases, such as hepatitis C and chickenpox, virus particles persist within hosts for many years. The set of values concentrated on in this paper were based on our earlier work, where non-oscillatory solutions of the differential equations were sought to avoid recurrence of disease by virtue of individual dynamics as opposed to interactive effects.

There are many potentially interesting extensions of the preliminary investigation of the model we have reported here. The path of the spread of infection could be followed graphically and with larger populations the existence of wave-like phenomena could thus be examined. One could examine the effects of randomly arriving disease-bearing immigrants on a population that has settled to an equilibrium situation with sub-threshold viral levels, described as a type (c) response at the beginning of Section 4. In addition, analysis could be performed on the effects of vaccination, quarantine and host deaths. Mutation of the virus to a more contagious form could also be studied by making the parameters p trans and the within host viral dynamical parameters take on new values for the members of the population who host the new mutant.

An examination of the Reed-Frost theory of epidemics

Infectious Diseases of Humans

World Health Organization, Epidemic and Pandemic Alert and Response

Transmission and dynamics of tuberculosis in generalized households

A multi-city epidemic model

A multi-species epidemic model with spatial dynamics

Limits of a multipatch SIS epidemic model

The Mathematical Theory of Infectious Diseases and its Application

Modelling the dynamics of LCMV infection in mice: conventional and exhaustive CTL responses

Mathematical model of antiviral immune response. III. Influenza A virus infection

Influenza: forecast for a pandemic

Mathematical Epidemiology of Infectious Diseases

Modelling disease outbreaks in realistic urban social networks

Smallpox as a biological weapon

The mathematics of infectious diseases

Stochastic models for epidemics: current issues and developments

Emergency response to a smallpox attack: the case for mass vaccination

Contributions to the mathematical theory of epidemics. Part I

Persistence of viral infections on the population level explained by an immunoepidemiological model

Variability in early HIV population dynamics

Transmission dynamics and control of severe acute respiratory syndrome

Spatiotemporal dynamics of epidemics: synchrony in metapopulation models

Strategy for distribution of influenza vaccine to high risk groups and children

The stochastic dance of early HIV infection

Hepatitis C viral dynamics in vivo and the antiviral efficacy of interferon-␣ therapy

Preparing for the next pandemic

Finding optimal vaccination strategies for pandemic influenza using genetic algorithms

Mathematical Models for the Growth of Human Populations

Preclinical assessment of HIV vaccines and microbicides by repeated low-dose virus challenges

Mathematical modeling of influenza virus replication in mammalian cells

Pandemic influenza preparedness: the critical role of the syringe

Semi-empirical power-law scaling of new infection rate to model epidemics with inhomogeneous mixing

Elementary Applications of Probability Theory: with an Introduction to Stochastic Differential Equations

Nature of equilibria and effects of drug treatments in some simple viral population dynamical models

Some properties of a simple stochastic epidemic model of SIR type

Spatial epidemic network models with viral dynamics

Enhancement of epidemic spread by noise and stochastic resonance in spatial network models with viral dynamics

Epidemic spread and bifurcation effects in two-dimensional network models with viral dynamics

An epidemic model in a patchy environment

Population-wide benefits of routine vaccination of children against influenza

The importance of lytic and nonlytic immune responses in viral infections

HCT thanks the Max Planck Institute for financial support and Prof. Dr. Juergen Jost for his kind hospitality. LT acknowledges financial support from A.C.I. Systèmes Complexes en Sciences Humaines et Sociales. We thank the four anonymous referees whose valuable criticisms and suggestions led to a much improved manuscript.