Estimating spatial coupling in epidemiological systems: a mechanistic approach
Go here for SFXEditor, M. Hochberg
Abstract
In recent years, ecologists and epidemiologists have paid increasing attention to the influence of spatial structure in shaping the dynamics and determining the persistence of populations. This is fundamentally affected by the concept of ‘coupling’—the flux of individuals moving between separate populations. In this paper, we contrast how coupling is typically implemented in epidemic models with more detailed approaches. Our aim is to link the popular phenomenological formulations with the results of mechanistic models. By concentrating on the behaviour of simple epidemic systems, we relate explicit movement patterns with observed levels of coupling, validating the standard formulation. The analysis is then extended to include a brief study of how the correlation between stochastic populations is affected by coupling, the underlying deterministic dynamics and the relative population sizes.
INTRODUCTION
Spatially structured models and emergent heterogeneity are becoming increasingly important features of ecological and epidemiological modelling ( 43; 36; 11; 23; 6). One of the most common forms of capturing the effects of space is the patch or metapopulation model (30; 2; 22; 37; 38; 21; 26). Such models consist of multiple subpopulations with movement (or ‘coupling’) between them, and hence provide an ideal framework for studying the dynamics of diseases within human communities (18).
Generally, increasing the strength of coupling between populations reduces spatial variation and acts to synchronize their dynamics. Sufficiently strong coupling may eventually lead to perfect synchrony or ‘phase‐locking’. The extent of this effect is, however, known to depend on the precise geometry of the coupling and the intrinsic patch dynamics (1; 3; 39; 7; 12). Spatial synchrony is a key ecological (and epidemiological) phenomenon since it can profoundly affect the likelihood of metapopulation persistence. In the absence of spatial synchrony we observe rescue effects, which prevent localized (subpopulation‐level) extinctions either from occurring or becoming permanent, thus increasing the probability of global persistence (10; 33; 13, 7; 19; 26).
Coupling between distinct populations can take a variety of forms. In most models of ecological population dynamics, coupling arises primarily from migratory dispersal (19). Hence, the important features are those determining individual movement behaviours, which are quantifiable (though the precise estimation of these is often very time‐consuming and rarely straightforward). In epidemiological systems, we are primarily concerned with the transfer of infection, which may occur via different routes. For sessile organisms such as plants, coupling is often due to external transport (e.g. by wind or insects) of the pathogen (24; 41). A similar mechanism of dispersal largely accounts for the spread of diseases such as foot‐and‐mouth disease or swine‐fever between farms (29). If one is interested in the epidemiology of childhood infections (such as measles or whooping cough), however, two important distinctions must be made.
Firstly, since these infections are directly transmitted, external transportation agents are not involved. Second, and more importantly, coupling no longer refers to the permanent translocation of individuals or infection from one centre to another. Instead, it is concerned with the relative levels of ‘mixing’ within and between populations, with temporary movement assuming a far more significant role. Unfortunately, good data on relevant human mobility patterns are hard to find. We do, however, have access to excellent data sets on the spatio‐temporal incidence of childhood infections (27; 35, 36). From these data we can estimate the correlation between epidemics in different populations, and attempt to infer possible movement patterns.
In this paper, we will use the classic SIR (Susceptible–Infectious–Recovered) framework for disease dynamics to investigate the role of movement between human populations. First we consider how models with simple phenomenological coupling can be related to more complex mechanistic models which explicitly describe the movement of individuals between populations. While this relationship takes a simple form when the two populations are of equal size, it is far more complicated given hierarchical populations. In the final section, we attempt to link models with data by relating variation in correlation with coupling strength and how this is affected by population size and the inherent population dynamics.
THE MODEL
We shall take as our basic framework the standard SIR model (4), where individuals are classified according to infection status (susceptible, infectious or recovered): 
where N (= S + I + R) is the total population size, b and d determine the per capita birth and death rates, respectively, and g−1
gives the infectious period. The parameter β represents the disease
transmission rate and is often assumed to be constant, giving rise to
damped oscillations in the deterministic model. In some of the best
studied cases (e.g. childhood microparasitic infections), there is
strong seasonal variation in contact rates and β is time‐dependent; this
may lead to more exotic dynamics (40; 16; 12; 28).
In most of the analysis that follows, we assume it to be constant,
although we comment on numerical findings for seasonally forced models.
We also set the birth and death rates to be equal (b=d ) such that the total population size, N, remains constant.
LINKING COUPLING AND MOVEMENT PATTERNS
In systems where a number of populations are coupled
together, the effect of spatial structure is commonly achieved by
allowing the ‘force of infection’ (the per capita likelihood of
contracting the infection; 4) in each population to be influenced by the level of infectious individuals in the other population (31; 8; 13; 36; 26). For example, for a system of two identical populations (N1=N2=N ), the spatially structured set of equations are given by:
The parameter σ defines a ‘coupling strength’ and is related in a phenomenological way to the movement rate between the two populations. Note that σ merely reflects how much ‘mixing’ exists between the two populations, measuring the proportion of infection that occurs between rather than within patches—it does not specify a dispersal rate. The precise formulation of Equation (2) is chosen such that the basic reproductive ratio, R0=β/g, remains constant as the strength of the coupling, σ, changes. This is a consistent feature of all the models given in this paper, and is vital if we are to reliably compare models with different couplings.
Such a method of coupling populations is clearly an over‐simplification of the true interactions between communities, although it has been widely used in the literature (17; 9; 31; 42; 36; 12). Implicit within this formulation is the idea that infectious individuals from population 1 can infect susceptibles in population 2, but this is achieved without any explicit mechanism for the transfer of infection. The primary aim of this paper, therefore, is to understand how basic movement patterns affect the overall ‘coupling’ between populations.
In order to construct a more mechanistic model, we need to specify the exact mixing patterns between the populations. The precise mechanisms by which infections are transferred to new centres are difficult to establish. Clearly, they involve the interchange of individuals between populations, but without precise data on patterns of human mobility, it is hard to describe what may be ‘realistic’. We shall initially consider the interaction between two populations of the same size. Individuals may leave their ‘permanent’ population and enter the other (‘temporary’) population, where they will remain for a given period of time. Within each population, we have SIR dynamics.
The nomenclature for this model involves identifying all individuals by two subscripts, such that Sxy refers to susceptibles from population x that are currently situated in population y. This allows us to explicitly model the movement of individuals back and forth between the two populations (Fig. 1)
and to quantify the effects of this movement on the epidemiological
coupling. In this mechanistic model, infection can only occur when both
susceptible and infectious individuals are present in the same
population. We assume infection operates as pseudo mass‐action (25; 32),
although if both populations are initially of equal size this
assumption has no bearing on the dynamics. The equations have three
separate components, births and deaths, infection, and movement between
populations; for example, the dynamics of Sxx are given by:
where ρ and τ are the rates at which individuals leave and return to their home location, respectively (as shown in Fig. 1).
The remaining equations are given in Appendix 1. For convenience, we
define and formulate results in terms of the parameter μ=ρ/τ, which
measures the ratio of time spent in the temporary to the permanent
population.
Schematic representation of the movement of infected individuals between two populations. Individuals in their permanent population (N11 and N12) visit the other, temporary population at rate ρ = μτ, those individuals in the temporary population (N12 and N21) return at rate τ.
In this (unforced) system, for most biologically realistic
parameters, the population levels rapidly converge to an equilibrium
such that the number of individuals in the permanent and temporary
location are related by:
Assuming this rapid population level convergence has already occurred, our mechanistic model (Appendix 1, equation 6) is an eight dimensional system, which we would like to equate with the much simpler four dimensional phenomenological model (equation 2). We achieve this by considering the eigenvalues of each system, and comparing how rapidly the populations converge in each formulation. In these unforced models, we always observe convergence of the two populations to a stable equilibrium point. By adding coupling between the populations the basic rate of convergence can be increased; and it is this increase that permits us to equate a level of coupling (σ) from equation (2) with mechanistic movement rates (ρ, τ ) from equation (6).
For all reasonable choice of parameters, the eigenvalues of the simple uncoupled and unforced model (equation 1) are complex conjugates, λ+0 and λ−0, with negative real parts. This corresponds to damped oscillations towards the fixed point (4). We expect the behaviour dictated by this pair of eigenvalues to dominate the behaviour of both the phenomenological and mechanistic models. Since the coupled system (equation 2) is four dimensional, it has two pairs of eigenvalues (λ+1,λ−1) and (λ+2,λ−2). One pair is always equal to the eigenvalues of the uncoupled system (without loss of generality, we assume that λ±2=λ±0) and provides information on the convergence to the fixed point. The second set of eigenvalues (λ±1) determine the speed with which trajectories in the two populations approach each other.
In Fig. 2(a), we plot the value of Real (λ±1−λ±0), which measures the extra rate of convergence of the two populations towards each other in addition to the normal convergence to the fixed point. The figure demonstrates that below some critical level (σc ≈ 0.1) greater coupling increases the synchronizing effect, and this is reflected by a more rapid convergence of trajectories in different populations towards each other. Above the critical coupling level, the eigenvalues λ±1 become real and so two values are observed in Fig. 2(a). In this region, increasing the coupling σ generally leads to slower convergence, as the larger eigenvalue dominates. Therefore, at the critical level of coupling, convergence and hence synchrony of the two populations is maximized.
(a) Comparison between the eigenvalues of the 4‐dimensional coupled system (equation 2), λ±1, and the eigenvalue of the 2‐dimensional SIR model (equation 1), λ±0. We note that Real(λ±1−λ±0) is always negative indicating faster convergence in the coupled model. For levels of coupling, σ, above 0.1 the eigenvalues λ±1 are real and hence two curves exist on the graph. (b) The relationship between occupancy ratio, μ, the relative return rate, 
and the phenomenological coupling, σ. Note that the coupling reaches a
maximum when μ=1 such that equal time is spent in both populations.
Asymptotic behaviour is reached whenever the movement rate τ is large
compared to the recovery rate g. (b=d=5.5 × 10−5 days−1, g−1=13 days, R0 = 
= 17).
The mechanistic model is described by a set of eight differential equations and therefore has four pairs of eigenvalues, which we denote by, Λ±n, (n=1, …, 4). Two of these pairs (say, Λ±3 and Λ±4) are large and negative; these correspond to the rapid distribution of individuals between the permanent and temporary populations, and can be ignored (for our purposes). Another pair of eigenvalues, say Λ±2, are equal to the eigenvalues of the uncoupled system (equation 1). This again leaves just one pair of eigenvalues remaining (Λ±1) to determine the effects of coupling on the disease dynamics. Comparing Λ±1(μ,τ) with Λ±1(σ), allows us to determine a relationship between the parameters of the two models.
Figure 2(b) shows the level of simple phenomenological coupling σ corresponding to the mechanistic parameters τ, g and μ. We note from the graph that the coupling is maximized when the returning rate (τ) is large (such that there is rapid mixing between the populations) and when individuals spend equal times in each population (μ=1). When τ is small, individuals spend a long time away from their permanent population, it is therefore quite likely that if they catch the disease they pass into the recovered class before they return home. Small values of τ mean that individuals cannot carry the disease between the two populations and hence there is a reduction in the amount of coupling.
It is interesting to note that in general as the infectious period changes, so does the relationship between μ, τ and σ. Thus, although the movement patterns are determined by the population and are independent of the disease, the epidemiological differences in the infectious period mean that the precise value of σ will be disease dependent. Analysis shows that it is the length of the infectious period relative to the average time spent away from home (τ/g) that determines the level of coupling (assuming that the infectious period is much shorter than the life expectancy). We notice that when the time spent in the temporary population is smaller than the infectious period τ−1 < g−1, the eigenvalues are very close to the large τ limit (see below). Intuitively, in most applied settings we might expect short stays in comparison to the infectious period so the precise value of τ and g are likely to be unimportant. (For infections like measles and whooping cough, we are dealing with infectious periods of 5 and 15 days, respectively, and an intuitive assumption would be that ‘mixing’ occurs on the scale of a day). Therefore, in such cases it may be appropriate to study the large τ limit, where precise disease parameters do not affect the level of coupling.
When the returning rate (τ) is very high, it is possible to
obtain a simple algebraic relationship between the strength of coupling
(σ) and distribution of individuals (μ). This is due to the rapid
convergence of the distribution of individuals between the two
populations such that,
Summing the differential equations for the number of susceptible individuals whose permanent home is population x gives:
with a similar equation for Ix. These are clearly identical to the phenomenological equation (2) with the strength of the coupling as given by
Hence in the majority of applied settings when τ is large compared to g, the full mechanistic model rapidly converges to the much simpler coupled model, and there is a simple relationship between the parameters of the two models.
A similar theory relating phenomenological coupling to movement parameters can be developed for the more realistic situation where the two coupled populations are of different sizes (Appendix 2). The analysis shows that given identical intrinsic movement rates in the two populations, coupling levels are identical except for a scaling factor determined by the relative population sizes.
CORRELATIONS FOR COUPLED STOCHASTIC MODELS
Although a detailed social survey may allow us to estimate the parameters μ and τ and therefore find a corresponding value of σ, for many diseases we are simply faced with the number of reported cases in each community. In this section, using a stochastic or Monte Carlo version of equation (2) (34), we compare the correlation between cases in two model populations, with the level of coupling σ. With this stochastic model, all events (birth, death, infection and recovery) occur at random, but with the same underlying rate as predicted by the differential equations. In addition, to prevent permanent stochastic extinctions of the disease, a small immigration rate (ɛ) of infectious individuals was added to the model.
For the time‐series of infectious cases from the two populations, I1(t) and I2(t) the correlation is defined as,
Here 〈·〉 corresponds to the long‐term average. Figure 3(a) shows the correlation (as calculated from a 5000‐year sample) for two populations of equal size N, with coupling σ.
The correlation between number of infectious
individuals in two populations using a stochastic (Monte Carlo) version
of the coupled model (equation 2). Part (a) clearly shows that the
correlation is largely unaffected by population size, and that the
simulation results are a close fit to the theoretical predictions (thick
lines). (b=d= 5.5 × 10−5 days−1, g−1=13 days, R0 = = 17, ξ = 0.0119, following the work of 5 the import rate was chosen as ɛ=5.5 × 10−5√Ν days−1).
Part (b) gives the value of ξ which provides the best fit to the
correlation, in terms of minimizing the mean squared error. The shaded
square corresponds to the epidemiological parameters used throughout the
rest of the paper. (b=d=5.5 × 10−5 days−1, σ=10−3 to 0.5, N=105.
Theoretical prediction
By formulating an equation for the covariance between the
number of infectious individuals in the two populations, we can begin
to understand the behaviour of the correlation. This is achieved using
moment closure approximations (34; 26)
in which means, variances and covariances are modelled but third order
cumulants are ignored. As detailed in Appendix 3, we obtain the
following approximation for the correlation,
where ξ is a parameter which can be estimated from numerical simulations. In Fig. 3(a)
, we find that ξ ≈ 0.0119, and is independent of population size—the
shaded curves show the theoretically predicted correlation using this
value. Thus, for a given disease and set of demographic parameters,
there is a unique curve relating the amount of correlation observed with
the strength of coupling. Hence, this result allows us to infer
individual level parameters from global behaviour. Given the disease
dynamics of two simple coupled populations, it is possible to calculate
the phenomenological level of coupling and therefore the individual
movement rates between the populations.
Figure 3(b) shows the value of ξ which provides the best fit to stochastic simulation results for a range of disease parameters. Again we find that equation (5) provides a good description of the correlation. The parameter, ξ, is large whenever both R0 and the infectious period ( g−1) are large; in this region of parameter space the correlation C increases slowly and almost linearly with the coupling. The increase in ξ for small values of R0 may be attributed to frequent extinctions within the model populations. In this regime the theoretical prediction is less accurate, but still gives a good qualitative description of the relationship between correlation and coupling.
Spatial complexity and seasonality
So far we have only considered the unforced system (which possesses a globally stable fixed point) and two isolated populations. In this section we shall consider the difficulties that arise in real applications when there is often strong seasonality and as well as coupling between numerous locations.
To capture the cyclic behaviour observed for many diseases, seasonality in the contact parameter β has to be included (14; 40). Here we use the formulation estimated for measles dynamics, where β takes a high value during term‐times and a lower value during school holidays (14; 15). This produces biennial dynamics in the deterministic model for measles parameters in developed countries (4; 28).
When seasonality is included, we need greater care in comparing the dynamics of the two populations. What we wish to measure is the correlation in the deviation from the average cyclic pattern. With measles parameters, for example, it becomes necessary to define a separate correlation Ct for every point t (0 ≤ t ≤ 2 years) on the biennial attractor. The correlation between two cyclic populations can then be taken as the average over the entire cycle of the separate correlations. With this modified definition of C, we find that the theoretical prediction (equation 5) remains a reliable approximation. We note that if the cyclic dynamics were ignored, and we simply calculated the correlation in the standard manner, then synchronization (or phase‐locking) of the epidemics in the two communities would produce much higher correlations than expected when the populations are large (N > 105 ).
This methodology must be used with caution. For some diseases, such as whooping cough (36) or rubella (4), seasonal changes drive irregular epidemics. Other diseases (for example influenza) are dominated by localized extinctions and subsequent recolonization, as well as other complications relating to antigenic variation. For such diseases, it is difficult to define an average cyclic pattern, and our theoretical approximation breaks‐down whenever deterministic epidemics become out of phase.
Most applied situations also require us to consider multiple populations with different degrees of coupling between them. A complete analysis of all the vast number of possible scenarios is not feasible. Here, we note in brief the results from two simple spatial configurations—arranging n identical populations in a line with either local or global coupling. For both these situations the same basic theoretical relationship between coupling and correlation still holds. With global coupling (every population coupled to every other population with strength σ) the value of ξ shows a very slight decrease with the number of populations. For local coupling the situation is more complex, ξ increases as we look at the correlation between ever more distant populations.
CONCLUSIONS
In human epidemiology, the interaction between discrete settlements can have a profound dynamic effect and important public health consequences. Most notably, re‐infection from an external source (rescue effects) plays a significant role in the global persistence of many diseases (17; 18). It is therefore important to have reliable models for the interaction between communities. We have established a clear equivalence between the standard phenomenological models of coupling, and a more mechanistic approach based on individual movement patterns. While this research has been primarily driven by patterns of human mobility, our approach should be applicable to any territorial organism, which spends short periods away from its permanent population.
The relationship between the phenomenological coupling and
the mechanistic movement rates can be calculated numerically for any
combination of disease and demographic parameters. It is interesting to
note that the level of coupling is sensitive to the length of the
infectious period, as well as on the movement rates. Despite identical
underlying movement behaviours, the coupling between communities is
disease dependent. This is because when the infectious period is short
and the length of stay in the temporary population long, it is unlikely
that an individual can acquire the disease in one population and return
home before recovering. However, for most human infections the time
spent out of the home population is short relative to the incubation
period of the disease. In this case, there exists an exact analytical
relationship between the phenomenological and the mechanistic models,
For most human populations, it is difficult to assess the
level of movement between different communities. Even when records of
travel do exist, they are rarely stratified sufficiently for us to
identify the potential mixing between susceptible and infectious
individuals. For childhood diseases, the distribution of disease is age
dependent, so information on the movement rates of recovered adults is
irrelevant; similarly for sexually transmitted diseases, the risk
associated with each individual is highly variable and not easily
ascertained. Often, however, we have extremely detailed reports of the
number of cases in each community. Stochastic simulations are used to
generate surrogate data for two coupled populations. Results from these
simulations agree with simple analytical predictions that in many
situations the correlation, C, is related to the coupling, σ, by
where ξ is a function of the particular disease parameters, but does not depend upon population size.
This work has mainly focused on the most tractable scenario, the interaction between two communities of equal size where the underlying disease dynamics possess a fixed point attractor. However, some consideration has also been given to the three main complications: different sized populations, cyclic dynamics and multiple interacting populations.
• When the populations are of different sizes, we again find that a phenomenological model is equivalent to the mechanistic form, although four distinct coupling parameters become necessary. If we again assume that visits to another population are of much shorter duration than the infectious period then a simplified analytical approximation with just two coupling parameters is achieved. Correlations between populations of vastly different sizes do not agree well with the theoretical formula, but it still provides a good approximation when the population sizes are within an order of magnitude of each other.
• If the dynamics have a regular cyclic pattern, then we must subtract this pattern from the time‐series before calculating the correlation; hence we are looking at the correlation between the deviations away from the attractor. The correlation then obeys the same formula as before, with a similar ξ value as a non‐cyclic (non‐seasonally forced) model. When diseases are strongly influenced by stochastic forces and epidemics are irregular and out of phase, it is difficult to assess the effects of any small degree of coupling.
• When many populations interact, the array of possible mixing patterns is vast. Concentrating only on local and global interactions, with equal coupling strengths, we firmly believe that the same relationship between correlation and coupling still exists, although the value of ξ may depend on the number of populations and the nature of the coupling.
Throughout we have assumed an SIR‐type formulation, so that individuals are infectious as soon as they are infected. However, for many diseases the SEIR model (4) is more appropriate; this introduces an exposed class when individuals are have caught the pathogen, but are incubating the disease and are not yet infectious. All of our results transfer in the obvious manner to this SEIR model, with very little change in the behaviour or quantitative predictions.
We have found clear analytical relationships between the individual movement patterns and phenomenological coupling, and then between coupling and the observed correlations. These relationships fairly robust, and therefore provide us with an important tool for analysing disease dynamics in a spatial context.
Acknowledgements
We thank the four anonymous referees for their comments on this paper. This research was supported by the Royal Society (MJK and PR).
Appendices
APPENDIX I
Mechanistic equations for populations of equal sizes
where x = 1 or 2 and y ≠ x. This formulation corresponds to individuals leaving their permanent population at rate ρ and returning at rate τ (Fig. 1).
Note that individuals inherit their parents’ permanent location,
irrespective of where they are born. Susceptible individuals can only
catch the infection from infectious individuals in the same location.
APPENDIX 2
Populations of different sizes
When the two populations are of different sizes then we necessarily have to consider that the parameter μ may be dependent upon Nx and Ny. However, it would seem sensible to assume that the time spent in the other community (τ−1)
is independent of the population sizes. To approximate the value of μ,
let us assume that individuals (in either population) randomly visit
other individuals (irrespective of location) at a rate γτ. The value of ρ
and therefore μ can be calculated from the rate of visiting an
individual in the other population. ρx is the rate at which someone in x visits someone in y, and is calculated as the visiting rate times the probability that the visit is to someone in y:
hence the value of μ is maximized when the other (temporary) population
is very large in comparison to the permanent population. We note that
the total number of individuals currently in a population Px=Nxx + Nyx,
is no longer equal to the number of individuals that belong to that
population; but again we expect rapid convergence of the distribution of
individuals between patches such that,
We can now reformulate the full model,
Such a formulation still retains the same value of R0
(irrespective of movement rates) and still allows us to reduce the
system to four dimensions, although the coupling terms are more complex.
In the limit where τ is large, the four coupling terms σxx, σxy, σyx and σyy are found to converge to two asymptotic values,
APPENDIX 3
Theoretical calculation of the correlation
Throughout this section we shall assume that the two
populations are of equal sizes, as this greatly simplifies the
mathematics and notation. Let CXY be the covariance between X and Y in the same population,
and C̄XY be the covariance when X and Y are in different populations,
Where 〈·〉 refers to the long‐term time average of a given quantity. We can now formulate an equation for the covariance, C̄II if we make a moment closure approximation and ignore the third order cumulants (34; 26).
If we make the simplifying assumption that IC̄SI/N is small, the differential equation has the fixed point,
From the equation for the number of infectious individuals, we find that
and therefore obtain the following equation for the correlation,
For simplicity we set,
and calculate the value of ξ from numerical simulations.
