Effect of Movement on the Early Phase of an Epidemic.

The early phase of an epidemic is characterized by a small number of infected individuals, implying that stochastic effects drive the dynamics of the disease. Mathematically, we define the stochastic phase as the time during which the number of infected individuals remains small and positive. A continuous-time Markov chain model of a simple two-patch epidemic is presented. An algorithm for formalizing what is meant by small is presented, and the effect of movement on the duration of the early stochastic phase of an epidemic is studied.

Introduction

Before infectious diseases become full-blown epidemics, they go through an initial phase during which they are subject to sometimes acute stochastic variations. From the perspective of the pathogen, this is a make-or-break phase, where outcomes could go either way: if transmission is to too few individuals, then the disease goes extinct, whereas if enough individuals are infected, an outbreak occurs.

Consider, for instance, the situation shown in Fig. 1, showing SARS-CoV-2 cases in Campbell County (Wyoming) in the early stage of the COVID-19 pandemic. We first observe short outbreaks, isolated in time, with up to 8 cases, then more and more frequent such outbreaks until the start of a major outbreak some time around 20 September 2020, after which case numbers undergo an exponential growth phase for some time (not shown).

Fig. 1

SARS-CoV-2 incidence in Campbell County, Wyoming, from the first case on 2020-03-20, to the start of the exponential increase in case counts (last date shown 2020-10-16)(Color figure online)

The first focus in this work is this pre-exponential phase; in particular, how long does it last and at what prevalence level does the switch occur between such behavior and exponential growth? Note that, as with most recent epidemic outbreaks, the intensity of the effort by public health authorities to limit disease spread must be taken into account: it is likely that the long period before the outbreak seen in Fig. 1 would have been shorter in the absence of public health interventions; see, e.g., Arino et al. (2020).

On the other hand, observe that in Fig. 1, until roughly 6 June, it appears that the disease could die out. It is likely, in this situation, that the observed persistence and transition into an epidemic was the consequence of case importations from nearby jurisdictions. This highlights the second question addressed in this work: does movement between jurisdictions have an effect on the duration of the stochastic phase?

The same type of initial stochastic phase can be observed in many datasets for a variety of different infectious diseases, see, e.g., Andersson et al. (2008), Brandt et al. (1973), Cowling et al. (2008), Hulth et al. (2009), Murphy et al. (1981) and Noyola and Mandeville (2008). One should also bear in mind that there are probably many such phases that do not lead to widespread epidemics and thus are not detected.

Whether these stochastic fluctuations are intrinsic or extrinsic is a hard problem. Altogether, it is indeed likely that the early stochastic phase of an epidemic is driven in part by a succession of “failed” epidemics (from the point of view of the pathogen), in which local extinction would take place were it not for a re-seeding by individuals importing the pathogen from different spatial locations.

There is also a mathematical basis for the distinction of the initial stochastic phase suggested by various approximation techniques developed in the literature. The deterministic and diffusion approximations of Kurtz (1970, 1971) are asymptotic in the population size P, but lose accuracy when the population or the population in any compartment is small; see Ovaskainen and Meerson (2010); Rebuli et al. (2017). This fact is exploited by the hybrid methods of Barbour (1975) and Rebuli et al. (2017). In the former case, branching processes are used in the region of the state space with a small number of infectives, while the CTMC is used in the latter. In both cases, the stochastic models are coupled with deterministic or diffusion approximation once the number of infectives becomes large. LATS (Milliken 2019) is another approximation technique, which constructs a finite state space with a small number of infectives around the disease-free quasi-stationary distribution to study first-passage problems including the probability of a minor epidemic. The algorithm employed below (see Appendix D) to choose the threshold number of infected individuals likewise seeks to minimize the error of calculating the probability of a minor epidemic.

In this paper, we focus on the initial stochastic phase evident from the data and previous mathematical results. We consider a two-patch metapopulation SIS epidemic model coupled by movement between patches, to investigate the effect of varying the overall rate of movement and the distribution of the population among the patches on features of the initial stochastic phase of an epidemic. In Sect. 2, we present the model and modeling assumptions and describe the methods used to formulate the stochastic phase and run the numerous simulations that form the basis of the figures presented throughout. The results of our numerical experiments are presented in Sect. 3. Our results are discussed in Sect. 4. In addition, standard analysis and results along with additional support for the choice of the threshold of infected individuals is presented in the appendices.

Model and Methods

In this section, we outline the model, its parameters and underlying assumptions, and the methodologies used to study the effect of movement on the stochastic phase.

Model

The model comprises two local populations that are coupled together by the movement of individuals. In the absence of movement, each location becomes isolated from the other. In that case, the disease dynamics in each isolated location are described by an SIS epidemic model, which is a variation on a model in Allen and Burgin (2000). A flow diagram of the interactions included in the model is presented in Fig. 2. The model takes the form of a continuous-time Markov chain (CTMC), details of which can be found in Appendix A.

Fig. 2

Flow diagram of the two population SIS epidemic model. The per capita rates of flow between compartments are shown (Color figure online)

The model features patch-specific disease transmission rates and recovery rates , . Movement rates from patch X to patch Y are given by for and .

Methods

Movement Parameters

There are no demographic transitions in the model, so the total population is conserved. We let the total population be given by P. Informed by analysis of a limiting deterministic model (see Appendix B), we note that when we ignore the distinction between susceptible and infected individuals, there is a quasi-stationary distribution for the populations in each patch given by , where and depend on the movement rates between patches. To facilitate analysis, we define the parameter to be the proportionor, in other words, and . Solving for yields the following expression for in terms of k and In the numerical experiments presented in Sect. 3, we consider the movement parameters and k to investigate the effect of movement on the duration of the stochastic phase.

Other Parameters

In addition to the movement parameters, other compound parameters and parameter assumptions are important for our investigation. The basic reproduction number is defined as the average number of secondary infections generated, in its infectious lifetime, by a single infected individual in an otherwise totally susceptible population. Despite the fact that this is a statistic that depends on the related mean field system, the basic reproduction number, , has been shown to be a threshold value in stochastic systems as well as deterministic ones by Allen and van den Driessche (2013). The basic reproduction number for the two-patch model is given by,and its derivation is given in Appendix B.

In the absence of movement, we have already noted that the two patches decouple. The basic reproduction in the isolated patch X is given by:where , , and are the transmission rate, recovery rate, and total population in patch , respectively. These patch are not threshold parameters per se but are used to quantify the overall favorability for disease spread in each location.

We are particularly interested in the case of a medium-sized city and a smaller satellite town or city connected by the movement of individuals. We assume that the communities have relatively similar access to health systems, and we consider the case of a single infectious disease spreading among their residents. As a result, we assume thatwhich corresponds to an average recover time of 1 week. This recovery time is slightly longer than the infectious period for influenza A and slightly shorter than that of mild-to-moderate COVID-19 (Centers for Disease Control and Prevention 2022).

The Stochastic Phase

Stochastic effects dominate population dynamics when populations or subpopulations are small and give way to deterministic population dynamics as populations increase, as implied by the limit theorems of Kurtz (1970, 1971). The hybrid models of Barbour (1975), Rebuli et al. (2017) incorporate a stochastic phase when the population is small in at least one compartment. In the case of epidemic data, time series can be characterized as having a phase of seemingly deterministic growth leading to a major outbreak preceded and followed by periods of stochastic oscillation in the number of cases.

We consider the stochastic phase prior to the deterministic growth that leads to a major outbreak. We are concerned with the case that the number of infected individuals is small. We consider the stochastic phase to be over when either there are no remaining infected individuals, or when the total number of infected individuals is sufficiently large. The latter case is associated with transition to deterministic growth to an endemic equilibrium of an associated deterministic model. Therefore, sample paths that reach a sufficiently large number of total infected individuals are associated with a major outbreak, while sample paths that reach zero total infected individuals before are associated with disease extinction.

Accurately characterizing the stochastic phase is therefore equivalent to determining an appropriate choice for the threshold number of total infected individuals . Since reaching total infected individuals is associated with transition to deterministic growth, we assume that once the threshold is reached, a major outbreak occurs with probability 1. For an appropriately chosen , a CTMC with stopping criteria and would approximate the probability of a major outbreak with a high degree of accuracy. Without loss of generality, we instead consider the probability of extinction. The true probability of extinction can be estimated various ways. One method is to choose stopping criteria as described above, but let be equal to the sum of infected individuals at the endemic quasi-stationary distribution predicted by the associated deterministic model. Such a model can be simulated repeatedly, and the probability of extinction can be estimated by the proportion of sample paths that terminate in the extinction stopping criterion . However, this method comes at a high computational expense. Instead, we approximate the probability of extinction using multitype branching process (MTBP) approximation (see Milliken (2017); Tritch and Allen (2018)). Details of the MTBP approximation can be found in Appendix C. Since the duration of the stochastic phase is partly dependent on the choice of , we assumefor all parameter choices. The fact that in a population of one million individuals is consistent with the error bound provided by the LATS method of Milliken (2019). Additional supporting evidence and greater detail regarding the choice of can be found in Appendix D.

Numerical Simulations

The heatmaps in Sect. 3 are obtained by varying a pair of parameters while leaving the remaining parameters fixed, resulting in over 3000 parameter combinations. For each parameter combination, we simulate an ensemble of 100,000 sample paths. In some cases, sample paths wander into regions of the state space where the inter-event time becomes very small. This can be viewed as similar to stiffness in simulating ordinary differential equations and leads to high computational expense. To avoids these problems, we employ an adaptive tau leap algorithm for approximate stochastic simulation using R package adaptivetau. The algorithm is run in parallel on multicore research computing servers.

The Effect of Movement

As noted in Sect. 2.2.1, all movement is controlled by the rate of movement from patch A to patch B, and the long-term average proportion k of the total population present in patch A. Recall, (2) shows that and k uniquely determine . The duration of a minor epidemic has been studied by Tritch and Allen (2018). As a control, we first illustrate how the duration of the stochastic phase differs from the duration of a minor epidemic in a single population. We then investigate the effect of changing (respectively k) both alone and varying with other key system parameters.

Duration of the Stochastic Phase in a Single Population

The term minor epidemic appears in Whittle (1955) and is characterized as being generally shorter in duration and having fewer cases than a major one. In the context of what we call the stochastic phase, minor epidemics are sample paths that lead to the extinction stopping criterion without reaching the threshold . Tritch and Allen (2018) found that the mean duration of a minor epidemic peaks around the threshold criterion . However, in that paper, Tritch and Allen also note that when , sample paths that lead to a major outbreak have a longer mean duration, on average. In the stochastic phase, major epidemics are sample paths leading to the threshold before . Since the threshold is typically chosen to be smaller than the total number of cases at the endemic quasi-stationary distribution, the duration of these sample paths in the stochastic phase is typically shorter than their duration in an unabridged model. In Fig. 3, we take patch A to be a single patch isolated from patch B. We see that the mean duration of sample paths that lead to extinction (shown in green) peaks around , while the mean duration of sample paths stopping after reaching (shown in red) peaks for . As a consequence, the mean duration of the stochastic phase (which includes all paths and is shown in black) is shifted to the right of . The shaded regions show the range of durations across the simulations run for each value of . The light shaded region represents the sample paths in the 97.5 percentile and not in the 2.5 percentile. The dark shaded region represents sample paths in the 75 percentile and not in the 25 percentile (the interquartile range).

Fig. 3

Mean duration of the stochastic phase in a single population patch A for various values of (Color figure online)

The Effect of

Recall from (2), can be expressed as a function of and k. As a result,that is, we can vary the overall rate of moving between patches by varying . Thus, as approaches zero, movement between patches becomes a rare event and as becomes very large, movement becomes so frequent as to minimize local features and make the system resemble a homogeneous single population. Deterministic metapopulation models exhibit similar behavior as movement rates become arbitrarily large. Unlike deterministic models, however, our simulations show that for small, the behavior of the stochastic metapopulation model depends on the patch in which the disease is introduced. That is, stochastic metapopulations coupled by very rare movement between patches behave similarly to a system of isolated populations. Furthermore, as the strength of the coupling as measured by initially increases, the behavior of the model continues to depend on the patch of introduction. From the stand point of probability of extinction and mean duration of an outbreak, the model still behaves a little bit like a decoupled system, except that the features of each patch are perturbed by the influence of the neighboring patch. The behavior of a two-patch system as varies is presented in Fig. 4.

Fig. 4

Range of durations (shaded) and mean duration (line) of the stochastic phase versus with disease introduced in patch A (red) or B (blue) for , and (Color figure online)

Of course, is not the only parameter that effects the mean duration of the stochastic phase. A benefit of the simplicity of our model is that it is easier to investigate the correlation between the relatively few model parameters. The only model parameters besides and k are the patch dependent infection rates and the recovery rate . For each fixed k, (1) implies that and are also fixed. Since is fixed, varies linearly with , . However, since has a more natural interpretation, we represent varying with as varying with .

Figure 5 consists of heat maps of the mean duration of the stochastic phase varying and while keeping . The graph on the left shows the results for and the one on the right represents the case . The red lines are level curves for the mean duration of the stochastic phase, and the blue lines are level curves for . When , then . In that case, for each choice of , the mean duration is decreasing in . For each there is a value of which maximizes the mean duration and that value increases with . In the right-hand plot of Fig. 5, and , almost an order of magnitude larger. For small values of , the mean duration is decreasing in , at least for . This is likely due to the fact that so that, when is relatively small, infected individuals are making extended sojourns into a patch where the path to either extinction or the threshold is more direct. However, the mean duration’s dependence on is non-monotone. As increases in this scenario, the outsized increase in means that whenever an infectious individual moves from patch A to patch B it is much more likely to move back quickly. For large , movement becomes more frequent and individuals approach the deterministic behavior of spending approximately k proportion of their lifetime in patch A.

Fig. 5

Mean duration (in days) of the stochastic phase for varying values of and with where on the left and on the right. Blue: value of the reproduction number for the two-patch model as given by (3) (Color figure online)

Typically two-patch models exhibit symmetry and simplicity that make them less interesting to study. However, the introduction of the parameter k in order to preserve a roughly constant proportion of the population in each patch obscures this symmetry and reveals complex and interesting behavior in even this simple model. Figure 6 shows heat maps of the mean duration while varying and when and (left) or (right). Note that regardless of k, when is small, there is not much dependence on , since movement is rare and the disease is introduced into patch A. On the left, and the dependence of mean duration is on both and . This is due to the fact that is nearly an order of magnitude smaller than . On the right, where is nearly an order of magnitude larger, as increases, the dependence of the mean duration on is much stronger.

Fig. 6

Left: Mean duration (in days) of the stochastic phase for varying values of and with and . Right: Mean duration of the stochastic phase for varying values of and with and (Color figure online)

The Effect of k

In this section, we investigate the effect of changes to the population proportion parameter k on the mean duration of the stochastic phase. We have already seen some evidence of the effect of k in Figs. 5 and 6. In both those figures, there are two panels with on the left and on the right. Within each pair, we see that the simulation results are qualitatively similar for small values of . Differences between and begin to emerge as becomes larger. As previously noted, this is not directly the result of the population proportion itself, but rather its indirect role in determining .

Figure 7 also presents two heat maps of the mean duration of the stochastic phase varying k together on the left and on the right. For these illustrations, , , and the disease is introduced in patch A. When not varying, and . Both panels indicate that for (), the mean duration of the stochastic phase is almost independent of k. As k increases to 1, any movement away from patch A where the disease is introduced is followed with high likelihood by a quick return. That is, the model behaves even more like an uncoupled patch A system. It is evident from the panel on the right that, as k approaches 1, the mean duration of the stochastic phase rapidly shifts from strong dependence on to becoming nearly independent of the reproduction number in the neighboring patch.

Fig. 7

Mean duration (days) of the stochastic phase as k varies with (left) and (right)(Color figure online)

Comparing the Effect of and k

To this point, we have investigated the effect of varying and the effect of varying k on the mean duration of the stochastic phase, but we have yet to fully consider the effects of varying both movement parameters simultaneously. Still, the pairs of heat maps in Figs. 5, 6, and 7 do give us the basis to make the hypothesis that the mean duration is strongly dependent on changes to and only begins to depend on k when it is near enough to 1 so that becomes very large. To test this hypothesis, we prepared the heat map presented in Fig. 8. The mean duration is calculated for various values of and k with , , and . The figure confirms our hypothesis that the mean duration depends on and only on k near 1 when the value of k guarantees that becomes large relative to .

Fig. 8

Mean duration (days) of the stochastic phase as a function of and k (Color figure online)

Discussion

Data and mathematical theory agree that systems described by compartmental models exhibit stochastic behavior when the population of at least one compartment is small. However, as the population increases the behavior of such a system, and models trying to represent it, becomes more deterministic. In the case of epidemic models, stochastic effects are always a factor following the emergence or re-emergence of an infectious disease. The transient dynamics following (re)introduction of the disease before either extinction or the onset of deterministically driven dynamics is described as the stochastic phase. A formal method for characterizing the stochastic phase by defining a threshold number of infected individuals is proposed. The threshold number of cases, , is determined so that the resulting truncated model detects the probability of disease extinction to a desired sensitivity. It is interesting to note that in a system with a total population of one million individuals, is found to be a sufficient threshold number of cases across a wide range of parameters. Details regarding the choice of are presented in Appendix D.

The effect of movement on the mean duration of the stochastic phase is investigated numerically for the case of a relatively simple two-patch SIS epidemic model. The two-patch case was chosen to reflect the setting of a medium-sized city and a smaller satellite city with an infectious disease circulating between the two populations. Despite the simplicity of the model, numerical experiments still yielded complex and interesting results. The choice of only two population patches also served to highlight a key feature of stochastic metapopulation models that differs from deterministic metapopulation models taking the form of a system of nonlinear ordinary differential equations. Namely, dynamics are highly dependent on the patch of disease (re-)emergence.

The sample paths leading to the end of the stochastic phase may terminate in disease extinction or reach the threshold number of cases indicating an outbreak. On average, the duration of sample paths leading to outbreak is longer than those leading to extinction. As the favorability for disease spread, as measured by the reproduction number, increases, paths leading to outbreak become more likely, but all paths become more direct and duration decreases. Even in a single, isolated population, the mean duration of the stochastic phase peaks for to the right of the classic threshold . This is an important distinction, since the mean duration of paths that lead to extinction (also known as minor epidemics) peaks at .

The effect of movement is captured in parameters and k. The parameter represents the rate of movement from patch A to B, but also the overall frequency of movement. The parameter k represents the average proportion of the population in patch A, but also defines the rate of movement from patch B to A, for each value of . The effect of varying these parameters alone, together and with the patch-specific rate of disease transmission (represented by the patch-specific ), is studied by numerical experiment. For each parameter combination, simulations are performed and their duration averaged. The massive computational expense is mitigated by taking advantage of advanced stochastic simulation techniques provided by a tau leaping method in R (adaptivetau).

Analysis shows that the mean duration of the stochastic phase generally decreases with the rate of movement when the patch of disease (re-)emergence is coupled to a patch more favorable to the spread of the disease. This relationship is not monotone, however. For example, when k is large (greater than 0.5), increasing indicates an even larger increase in and the system resembles the patch of disease (re-) emergence being coupled to a patch less favorable to the spread of disease. In that case, the mean duration can actually begin to increase with increasing . In contrast, the mean duration of the stochastic phase does not seem to depend strongly on the proportionality constant k. However, this conclusion also has its exception when k approaches 1. The exception is driven by k’s role in defining . That is, even for fixed , increases to infinity as k increases to 1.

Analysis and simulations presented here are not meant to be interpreted as a real-world case study. Rather, they are meant to provide understanding of general changes to the qualitative behavior of the early stochastic phase of an epidemic resulting from changes to control the rate of disease spread locally or to control movement between populations. The stochastic phase is important to understand because it is the phase during which any particular outbreak could go away or lead to an epidemic. While some individual sample paths spent more than a year in the stochastic phase, the mean duration is generally less than a month. This means that public health strategies have to be able to be quickly deployed in order to have an effect on an outbreak during this critical phase. Additionally, the threshold indicates the need for ongoing development of rapid and reliable surveillance programs to catch an outbreak before it moves into a deterministic growth phase.

Table 1

Transition rates for the CTMC

Description	Transition	Rate
Infection in A	S^A,I^A,S^B,I^B→S^A-1,I^A+1,S^B,I^B	β^AS^AI^A
Recovery in A	S^A,I^A,S^B,I^B→S^A+1,I^A-1,S^B,I^B	γ^AI^A
Movement of S A to B	S^A,I^A,S^B,I^B→S^A-1,I^A,S^B+1,I^B	m^ABS^A
Movement of I A to B	S^A,I^A,S^B,I^B→S^A,I^A-1,S^B,I^B+1	m^ABI^A
Infection in B	S^A,I^A,S^B,I^B→S^A,I^A,S^B-1,I^B+1	β^BS^BI^B
Recovery in B	S^A,I^A,S^B,I^B→S^A,I^A,S^B+1,I^B-1	γ^BI^B
Movement of S B to A	S^A,I^A,S^B,I^B→S^A+1,I^A,S^B-1,I^B	m^BAS^B
Movement of I B to A	S^A,I^A,S^B,I^B→S^A,I^A+1,S^B,I^B-1	m^BAI^B