Dynamical response of multi-patch, flux-based models to the input of infected people: Epidemic response to initiated events.

The time course of an epidemic can be modeled using the differential equations that describe the spread of disease and by dividing people into "patches" of different sizes with the migration of people between these patches. We used these multi-patch, flux-based models to determine how the time course of infected and susceptible populations depends on the disease parameters, the geometry of the migrations between the patches, and the addition of infected people into a patch. We found that there are significantly longer lived transients and additional "ancillary" epidemics when the reproductive rate R is closer to 1, as would be typical of SARS (Severe Acute Respiratory Syndrome) and bird flu, than when R is closer to 10, as would be typical of measles. In addition we show, both analytical and numerical, how the time delay between the injection of infected people into a patch and the corresponding initial epidemic that it produces depends on R.

Introduction

Understanding the spread of infectious diseases through a population is important in determining the risks and consequences of natural or induced epidemics, such as those that are a consequence of infected people moving into a new area. An important challenge for these models is to incorporate the spatial heterogeneity in the geographical distribution of people in order to provide a more realistic account of the spread of the infection.

There are different ways to model the spatial-temporal patterns in a continuously distributed population. One approach is to use the set of nonlinear parabolic partial differential equations (PDEs) which incorporates both temporal and spatial processes at the same time [1], [2], [3]. PDEs are used to model a variety of ecological phenomenon such as the dispersal of species, ecological invasions, the effect of habitat geometry, and the formation of diffusion-driven spatial patterning [1]. These methods have shown that the movement of organisms can produce large-scale patterns in homogeneous environments and that the movement of multiple species can change the outcome of competition or predation in heterogeneous environments. Recently, the PDE approach has been used to model a spatial epidemic following the point release of a rapidly dispersing infectious agent. It was demonstrated that the resulting epidemic exhibits two distinct phases: the primary phase where the epidemic wavefront propates at constant speed and a secondary phase with a deceleration wavefront [3]. This approach can also be applied to understand the traveling waves of epidemics, such as those that have been observed in the spread of Dengue haemorrhagic fever [4]. Another approach is to use agent based models that follow each individual in the population [5], [6], [7]. Agent based models are either highly computationally demanding because they must include the several parameters of each individual agent, which for some models can consist of millions of agents, or they must assume some distribution of parameters across the set of agents without full empirical knowledge of the parameters of each agent.

Here we take a different coarse grained approach that captures important aspects of the spatial heterogeneity, yet is computationally simple, and can model specific cases of the spatial-temporal spread of disease. This approach is based on dividing the population into individual locations, called “patches”. The spread of infectious diseases in each patch can be described by ordinary differential equations [8], [9]. In the simplest classical model, we model the number of susceptible people, S, and infected people, I, homogeneously distributed in each patch, by using the following set of ordinary differential equations where N is the total number of people, μ is the birthrate at which new susceptibles are added to the population, is the contact rate, and γ is the recovery rate. Since we are concentrating on the transient response rather than long time behavior, these two equations are adequate and we do not need to include the death rate. The term in Eq. (1) is called the mass action term. Here we use the form of this term with a dependence, sometimes described as the “frequency dependent transmission” mass action term. Different dependencies on N for this term correspond to different assumptions as to whether diseases are spread proportionally to the number of infected people, the fraction of infection people, or some scaling that depends on the size of the population itself and how well mixed are the susceptible and infected people in that population. It is not clear what dependency on N is the most realistic one for this term. The form originally proposed by Hamer [10] and Kermack and McKendrick [11] was βIS which was based on an analogy to the kinetic rate constant of a first order chemical reaction between two well mixed chemical species. However, even chemical reactions, which become less well mixed as the reaction proceeds, display a different (and time dependent) dependence on N [12]. The review by McCallum et al. [13] concludes that “increasingly, the weight of evidence is that simple mass action is not an adequate model in many situations. A clear default alternative has yet to emerge”. With this uncertainty in mind, we have here chosen to use the form of this term, which we have found useful in our previous studies [14], [15].

This simplest form of the classical model having only one patch is not adequate to represent the spatial heterogeneity in the distribution of susceptible and infected people that is actually found. We therefore introduce the spatial heterogeneity by modeling the population as organized into a number of separate patches of different sizes with different infectious parameters in each patch [5], [6], [7]. We then compute the epidemics in each patch as the patches interact with each other. This patch approach provides a tractable way of modeling and studying the coarse grained spatial heterogeneity. We can then use these patch models to understand how the epidemics in a single patch are driven by its interaction with other patches and the injection of infected people into one of the patches.

We are particularly interested in understanding how the migration of infected people into a patch, or the movement of infected people from one patch to another, effects the spread of disease as these results shed light on the epidemics developed in response to natural occurrence or a purposeful initiated event.

In previous multi-patch models the interaction between the patches was described by a second order bilinear term formulated as an extension of the mass action term for the homogeneous population in one patch. It is useful to treat the interaction between the patches as a flux term in a more realistically, namely by modeling the direct movement of susceptible and infected people between the patches. Liebovitch and Schwartz [14] derived the interaction between the patches in this form and showed that this formulation could also lead to useful new analytical and numerical results. Here we extend their study of multi-patch, flux-based models by including: (1) the annual driving term for the infection which has been shown to be very important in describing the nonlinear properties of epidemics [15], (2) the injection of infected people into a patch, and (3) more complex geometries of connections between the patches.

In particular, we study how the epidemics induced by the annual sinusoidal driving term and the injection of infected people into one patch are transmitted to the other patches to which it is connected. These results have relevance to understanding the time course in response to a natural or purposeful initiated event. We also show how this dynamical response is different for diseases, such as SARS (Severe Acute Respiratory Syndrome) or bird flu which have a low reproductive rate compared to diseases, such as measles, which have a high reproductive rate. In addition, we find the delay time between an injection of infected people into a patch and the resulting first outbreak of an epidemic depends on the reproductive rate of the disease.

Based on the approach of Liebovitch and Schwartz [14] we use Eqs. (3), (4) with the variables and , which are where the flux of people during a time interval Δt from patch k to another patch j is given by . Since the fraction of infected people is quite low when there is no epidemic we introduce the logarithmic transformed variables so that they have values that are not close to zero when there is no epidemic. This helps increase the numerical accuracy of the numerical integrations. With the logarithmically transformed variables and , these equations now become where the flux of people during a time interval Δt that move from patch k to another patch j is given by and the number of individuals in each patch, remain constant.

The strength of the annual driving term has a strong effect on the dynamical behavior [15]. For this term we use the form where is the strength of the annual driving term and ω is the frequency. The can induce either periodic or chaotic dynamics. We investigated both cases. In order to determine an appropriate set of parameters for these models, we first computed the bifurcation diagram of independent patches for δ between 0.0 and 0.2. We then chose two representative regimes: (i) periodic where and (ii) chaotic where . We also set the other parameters as: , , , or 1575 yr−1, , r (fractional migration rate) in the range from 0.01 to 0.5, (the time at which the infected people are injected into one patch) = 50 yr and (birthrate of new susceptibles) = (birthrate of new susceptibles) = 0.02 yr−1. The reproductive rate defined as the number of new infected cases that are generated by each infected person, is given by , which for the values of or 157.5 yr−1 and yields and 1.575.

Two patches, intrinsic epidemics

In order to understand the dynamics of epidemics in an isolated large patch when infected people are injected into the patch, we first used a model with two patches, where repeated epidemics in a small patch injects infected people into the large patch. We studied two different dynamical regimes, where the intrinsic dynamics of the large patch was either periodic or chaotic. We used the parameters: , , , , , , and . Fig. 1 shows the fraction of the susceptible and infected people in the large patch (solid line) and in the small patch (dashed line). The unperturbed base dynamical state is shown up to the time during which time there is no flux of infected people between the patches . At the flux of infected people starts from the small to the large patch . Four cases are illustrated in Fig. 1 for the periodic and chaotic regimes each with and . We first consider the periodic regime for both values of R. Fig. 1(b) shows that for , the injection of infected people from the small patch induces corresponding epidemics in the large patch. After time, that we call the “infected-induced transient time”, the system to returns to its base dynamical state before the infectious people were injected. Fig. 1(c) shows that for , the injection of infected people from the small patch also induces corresponding epidemics in the large patch. However, unexpectedly, there are also additional epidemics, which we call “ancillary” epidemics, in the large patch even though the small patch has no epidemic at that time. These ancillary epidemics are not the result of new infected people entering the larger patch. What is happening here is that the first epidemics in the large patch pushes that patch into a new part of phase space, far from its original base dynamical state, causing the large patch to produce a second epidemic all by itself. We also studied how the presence of these ancillary epidemics depends on the rate of migration, r, from the small patch into the large patch. As r is decreased from 0.5 to 0.01, these additional epidemics decreased and the case then behaved as the case.

Fig. 1

Fraction of susceptible and infected people in each of two patches computed from the two patch, flux-based model when the flux of infected people from the small patch to large patch is initiated at in (b), (c) (d) and (e). The 4 cases shown are for the periodic regime with the reproductive rate (a) R=15.75 and (b) R=1.575 and for the chaotic regime with the reproductive rate (c) R=15.75 and (d) R=1.575. The flux rate is , . The solid line shows the fractions in the large patch and the dashed line the fractions in the small patch. Squares mark the largest number of infected people and susceptible in each patch.

In order to confirm that these additional epidemics in the large patch were caused by its perturbation from its base dynamical state we investigated what happens when an isolated patch is started with initial conditions off their base state values by a fraction between 0.001 and 0.9. When that isolated patch is started with the initial conditions between 0.1 and 0.68 off its base dynamical state values, it produces ancillary epidemics similar to the case of . This supports our conclusion that these ancillary epidemics in the case of are due to the fact that the epidemics from the small patch have sufficiently disturbed the large patch from its base dynamical values to produce these ancillary epidemics.

The chaotic regime is illustrated in Figs. 1(d) and 1(e). Fig. 1(d) shows that for , there are naturally occurring epidemics before the migration between the patches begins at . The injection of infected people from the epidemics in the small patch also triggers additional epidemics in the large patch, but these additional epidemics are relatively small compared to the naturally occurring epidemics in the large patch. Fig. 1(e) shows that for , the large patch has not only its naturally occurring epidemics, but also small epidemics driven by the smaller patch. These additional epidemics in the large patch are considerably larger than those induced in the case of .

In the subsequent figures we concentrate on the periodic regime, since the chaotic regime always produces naturally occurring epidemics which mask the epidemics induced by the injection of infected people. However, we also computed all these models for the chaotic regime and we refer briefly to those results in the text.

Infected input into isolated patches

We now extended the models to determine what happens when infected people from the outside are injected into a patch. To do this we add an additional term to Eq. (6) of the form where the variance is 1 yr and the mean number of infected people that are injected over that time, , is between 10 and 104. This form does not imply stochasticity, but was chosen as a simple form to describe a broad pulse of infected people introduced into the patch.

We first considered what happens when infected people are injected into isolated patches. We actually did this by using the two patch model with no migration of infected people between the patches. We consider the case of two unconnected patches in the periodic regime with either or . We varied amount of the number of injected infection people, , from 10 to 104. The population, , of large patch was 106. For the specific case shown in Fig. 2 , was 103. The other parameters are: , , , , , and .

Fig. 2

Infected fraction of people when infected people are injected at into two isolated patches in the periodic regime with the reproductive rate in (b) R=15.75 and (c) R=1.575. The solid line shows the infected fraction in the large patch and the dashed line the infected fraction in the small patch.

The results computed from this model are shown in Fig. 2. Fig. 2(b) shows that for , the periodic regime is well established before the injection of infected people injected at . This input produces a single large epidemic and then the patch returns quickly to its previous periodic behavior. But in Fig. 2(c) for , repeated epidemics ancillary are generated. As time goes by, the amplitude of these epidemics decreases, but this transient behavior lasts a considerable time before the patch returns to its previous base dynamical state behavior.

We also investigated the chaotic regime for and . In both cases there is one peak due to the input of infected people and the rest of the behavior is analogous to that found and illustrated in Figs. 1(d) and 1(e).

Infected input into two connected patches

Next we studied the response of a connected set of patches to the injection of infected people. The first case was a two patch model in the periodic regime where infected people are injected at into the small patch and there is a migration of infected people from the small to the large patch as illustrated in Fig. 3 (a) (, ). Fig. 3(b) shows that for , the injection of infected people into the small patch produces a transient increase (as expected) in the infected fraction in the small patch (dashed line), and a (unexpected) decrease in the fraction of the infected population in the large patch (solid line), and that both patches return quickly to their previous periodic behavior. On the other hand, Fig. 3(c) shows that for , the injection of infected people into the small patch generates epidemics both in the small and the large patch. There are also, as found in Fig. 1(c), additional ancillary epidemics in the large patch when there is no corresponding epidemic in the small patch. The transient behavior in both the small and large patches continues much longer for the case shown in Fig. 3(c) than for the case shown in Fig. 3(b).

Fig. 3

Infected fraction of people for two patches in the periodic regime when infected people are injected at into the small patch with the migration of infected people from that small patch to the large patch with reproductive rate (b) R=15.75 and (c) R=1.575. The solid line shows the infected fraction in the large patch and the dashed line the infected fraction in the small patch.

Infected input into radially connected patches

We then studied models with a large central patch connected to other satellite patches, roughly corresponding to a central urban hub and its surrounding suburbs. Such a model, with 6 patches connected to a central hub, in the periodic regime, is shown in Fig. 4 (a). The infected people were injected into only the largest patch (solid line) at . The migration rate, r, between the largest patch and the other patches was 0.1 and was bidirectional. Fig. 4(b) shows that for there are simultaneous epidemics in each patch which then return quickly to their previous periodic behavior. On the other hand, Fig. 4(c) shows that for there are slightly longer delays in the epidemics amongst the patches and that repeated epidemics continue for a long time. In this case, the infected-induced transient time is very long before the patches return to their original periodic behavior. In both of these models the infected people were injected into the largest patch. We also studied the result of infected people injected into the smaller satellite patches. In that case, as the satellite patches have much fewer people than the central hub, their effect on the central hub is quite limited.

Fig. 4

Infected fraction of people for a model with a central hub of one large patch connected to 6 smaller satellite patches with bidirectional migration of infected people between the hub and the satellite patches. All the patches are in the periodic regime. Each patch has 1/5 the number of people as that of the next largest patch. At infected people are injected into the large central patch. The reproductive rate is (b) R=15.75 and (c) R=1.575. The black line shows the infected fraction in the largest central patch.

Infected input into serial connected patches

We also investigated what happens when the patches in the periodic regime are serially connected and the infected people were injected into the largest patch as shown in Fig. 5 (a). First, we consider the cases with . In the model shown in Fig. 5(b), with 10 patches and , the infected populations in all the patches have epidemics which have different magnitudes in each patch. There is only a brief phase delay in the timing of the epidemics in the series of patches. This is because the epidemics along the line of patches are not triggered by epidemics in each patch which are then passed on to the next patch. Rather, it is the small influx of infected people which moves rapidly through all the patches which excites each patch, almost simultaneously, off of its base dynamical state and into an epidemic. In the model shown in Fig. 5(d) with 3 patches and , all the patches have one large simultaneous epidemic triggered by the injection of infected people. In both models, the infected-induced transient time is rapid and the system returns to its base periodic dynamical state is a short time. In the model shown in Fig. 5(c) with 3 patches , as in all the other cases when , the infected-induced transient behavior continues for an extended time before the patches return to their periodic behavior. The behavior of models with 10 patches and (not shown) was also similar to that with 3 patches and . We also studied the result of infected people injected into smaller patches in the series of patches. As in the case of the large central hub and satellite patches described above, those smaller patches have much fewer people than the large patch and so their effect on the large patch is quite limited.

Fig. 5

Infected fraction of people for models with 3 to 10 serially connected patches with bidirectional migration of infected people between adjacent patches. All the patches are in the periodic regime. At infected people are injected into the largest patch. Each line in the graph corresponds to different patch. In (b) there are 10 patches with R=15.75, and each patch has 1/10 the population of the next largest patch; in (c) there are 3 patches with R=1.575, and each patch has 1/5 the population of the next largest patch; and in (d) there are 3 patches with R=15.75, and each patch has 1/10 the population of the next largest patch.

Dependence of the time to the first epidemic on the reproductive rate

A salient feature from all these results is that the both the infected-induce transient for the system to return to its base dynamical state behavior, and the time to the first epidemic increases as the reproductive rate R approaches 1. We used both numerical simulations and analytical approximations to better understand this result concentrating on determining the time to the first epidemic.

Numerically, we studied the duration of the time to the first epidemic in models of 3 serially connected patches in the periodic regime. We estimated the duration of this time as the interval between the time at when infected people are injected into the largest patch and the time when the first corresponding epidemic is produced. Samples of these calculations for , 4.575, and 6.575 are shown in Figs. 6(a), 6(b), and 6(c) and the time to first epidemic outbreak as a function of R is shown in Fig. 6(d).

Fig. 6

The duration of the time to the first epidemic depends on the reproductive rate R. Shown here rate the infected fraction of people for models with 3 serially connected patches with bidirectional migration of infected people between adjacent patches and each patch has 1/10 the population of the next largest patch. All the patches are in the periodic regime. At infected people are injected into the largest patch. As can be seen with (a) R=1.575, (b) R=4.575, and (c) R=6.575, the duration of the time to the first epidemic after the injection of infected people decreases as R increases. (d) The functional form of the duration of this time as a function of R measured from the numerical simulations (Numerical, solid line) is a good match to that estimated from the analytical approximation in Eq. (12) (Analytical, dashed line).

Analytically, we can also estimate the duration of this time by studying the dynamical behavior of a single patch. Its trajectory has two different time scales: a slow time scale corresponding to the almost linear slow growth of the number of susceptible people and a fast time scale corresponding to the exponential growth in the number of infected people at the onset of the epidemic.

For the case of the slow time scale in the susceptible population, while , Eq. (1) can be approximated as The transient time to the epidemic at time T is the interval in time between 0 to T. Integrating Eq. (8) we get where the susceptible people in time series are increasing linearly while the infected people are in stable state with . The onset of the epidemic is determined from Eq. (2), where The number of infected people will grow exponentially when the number of susceptible people S has increased so that λ increases from being less than zero to being greater than zero, which from Eq. (11) implies that Since the onset of the epidemic at , Eq. (9) implies that Therefore combining Eqs. (12), (13), we find that since . This result is illustrated in Fig. 6(d) which compares too very different mathematical approaches, each with somewhat different limitations in their accuracy; namely, numerical finite difference integration and a global analytical approximation. Considering the variability of epidemics and the approximations in measuring them in the numerical computations and the approximations made in the analytical derivation, Eqs. (8), (9), (10), (11), (12), (13), (14), we note that both the analytical and numerical results have very similar functional forms, even though they are not equal. We consider the correspondence between the numerical and analytical form of this relationship valuable confirmation of the essential concept, namely, that the time to the first epidemic depends significantly and inversely on their reproductive rate R.

Summary

The multi-patch, flux-based equations provide a simple and effective way to study the base dynamical state and transient behavior of the spatial-temporal spread of disease in populations that are divided into different regions with the movement of infected people between those regions. We showed here how the intrinsic epidemics driven by an annual driving term and the injection of infected people into a patch can spread to other patches. The strength and timing of these subsequent epidemics depends on the strength and geometry of the migration between the patches and on the reproductive rate, R, of the disease. The most salient observation of the different conditions described here is that when R is close to 1, as is the case for SARS and bird flu, then the transients generated in the patches are of long duration and there are additional “ancillary” epidemics that are not due to further input of infected people, but rather are a manifestation of the fact that the previous epidemic has pushed the dynamical state of the patch far from its base dynamical state. These ancillary epidemics may have important implications for policy makers deciding how to respond to additional epidemics after an initial induced event. In this regime of and periodic parameters, the initial epidemic is produced by injection of infected people. However, the subsequent ancillary epidemics do not correspond to the injection of additional infected people. That is, the ancillary epidemic represents a rebound of the system, rather than a second event of injected infected people. In addition we showed both analytically and numerically how the time to the first epidemic after the injection of infected people depends on the reproductive factor R.