Calculation of Epidemic Arrival Time Distributions using Branching processes

It was said by Jules Verne in 1873 that “The world has grown smaller, since a man can now go round it ten times more quickly than a hundred years ago,” that one could travel around the world in eighty days [Verne, 1873]. One can’t help but wonder what Verne would think of the advances in technology that have occurred in the century and a half since, which now allow the circumnavigation of the globe in less than eighty hours. The World Aviation Network (WAN) has, over the past century, shrunk the globe to a fraction of its former size, allowing for vast increases in tourism, immigration, and trade. At the same time, diseases that might once have travelled at the speed of cart, ship or train now cross the world in a number of days or hours. In 2003 Severe Acute Respiratory Syndrome (SARS) was able to spread from China to Vietnam, Hong Kong and Canada within two weeks of being reported to the World Health Organization [WHO, 2003]. Similarly, the 2009 H1N1 flu pandemic first reported in Mexico was able to reach both Europe and Asia within a fortnight [WHO, 2009]. Now, in 2020 we see that coronavirus SARS-CoV-2 has spread quickly across the world [WHO, 2020], reaching all but a handful of countries withing a few short months.

Modelling of disease spread and disease arrival time is critical, both to mitigate the effects of a disease and to determine where best to apply quarantine and screening protocols so as to minimize the spread of infection. In order to model the spread of disease, given the current WAN, a number of highly detailed epidemic simulation models have been created, such as GLEAMviz[Broeck et al., 2011]. Depending on how they are set up, these models can include such factors as vaccination, incubation times, multiple susceptibility classes, non-linear responses, seasonal forcing, quarantine, and the stochastic movement of individual agents, providing a detailed and comprehensive modelling framework. Unfortunately, in many cases, such models quickly become black boxes in and of themselves - objects to be studied and analyzed only via costly simulation - allowing predictions of the future that can be observed, but seldom understood. This difficulty is further exacerbated by the difficulties of determining parameter values, often a delicate task during the initial stages of a disease outbreak when information is scarce.

Depending on the questions being asked, explicit modelling of such intricate details is necessary, however, for other questions, such as the determination of epidemic arrival time (‘AT’), it seems that simpler methods may suffice. For this reason, a number of authors [Brockmann and Helbing, 2013, Iannelli et al., 2017, Gautreau et al., 2007, Chen et al., 2018] have proposed a variety of heuristics and metrics – artificial measures of distance based on flight data from the WAN. The goal of such metrics is to predict relative arrival times for an epidemic starting in a specified location, predicting for example that an epidemic starting in Vancouver will take twice as long to reach Istanbul as it will to reach Rio de Janeiro. Absolute AT will also depend on the parameters of the particular disease in question.

A particularly elegant example of such a metric is Brockmann & Helbing’s [Brockmann and Helbing, 2013] ‘effective distance’, in which they define the direct distance between a pair of locations as

where $P_{ij}$ is the probability that a particular individual who is leaving location $i$ is travelling to location $j$. The total effective distance between nodes $a$ and $b$ is then given as the shortest path between a given pair of nodes:

where here we minimize over all possible paths from $a$ to $b$, and $d_{ij}$ denotes the effective distance along each step of this shortest path (and $P_{ij}$ the corresponding probability).

Papers by Gautreau, Ianelli, et al. provide more detailed and mathematically justified metrics, by first calculating the expected time for a growing epidemic to make its first ‘jump’ [Gautreau et al., 2007], and then summing over all possible paths joining a pair of airports (in an appropriate manner) [Iannelli et al., 2017]. More recent work by Chen et al. [Chen et al., 2018] made use of linear spreading theory and the matrix exponential, in order to attain similar accuracy at reduced computational cost.

Underlying each of the above efforts, either implicitly or explicitly, is an assumption of unbound growth; the idea that in the early stages of an epidemic, disease spread is not limited by population size. Such unbound growth can, and often has been, modeled as smooth exponential growth, a continuous approximation of discrete real world populations. An alternative model for unbound growth, one often used when small populations, rare events, and probabilities are of interest, is the branching process [Kimmel and Axelrod, 2015]. The branching process is a classical tool in the study of extinction and evolution[Schaffer, 1970], and has been used to study epidemics[González et al., 2010], surnames and genealogies [Watson and Galton, 1875], and cancers [Goldie and Coldman, 1979, Coldman and Goldie, 1986]. In this paper we apply the branching process framework in order to calculate not only the mean, but also the full distribution of ATs for arbitrary networks.

In section 2 we introduce a basic ‘multi-compartment branching process’ model, and from this model derive a system of ODE’s that precisely determine the probability of epidemic arrival by a given time. In section 3 we determine the mean and variance in the AT, and explore how these results both support and contrast with the results of past papers. Section 4 explores the sensitivity of our predictions to perturbations to system parameters and network structure. In section 5 we use data from real world flight networks, and compare our predictions, along with the predictions of past authors to epidemic arrival times as observed for in the 2003 SARS epidemic, and 2009 H1N1 influenza pandemic. In both cases we observe correlation between predicted and observed results, but also significant noise. We are, unfortunately, unable to reproduce previously published results, and instead observe the observed stochasticity is both larger than previously reported, and larger than might be expected given the underlying model.

Our goal in this paper is to provide a ‘canonical’ approach to calculating arrival times; that is to say, an approach which relies on minimal assumptions and approximations, and can be trusted to provide accurate results across all of parameter space. To the extent that the method works, it can be seen as providing a foundation to past methods and heuristics. To the extent that our predictions disagree with real world data this is suggestive of either flaws in the data, or gaps in the underlying model, gaps that will require not simply better mathematics, but instead better understanding of the system under study.

Let us begin by concretely defining the model; we consider a network of $N$ connected nodes, each node representing a location. In our case these nodes refer to particular airports in the WAN. Individuals travel from node $i$ to node $j$ at transport rate $\gamma T_{ij}$. Here $T_{ij}$ gives the relative transport rate along a given route, while $\gamma$ is the global flight rate. While $\gamma$ and $T$ can be folded together as one parameter, maintaining this distinction will prove convenient later. We select $T_{ii}$ to represent the rate at which individuals leave location $i$; $T_{ii}=-\sum_{j\neq i}T_{ij}$. Each node contains a population of $I_{k}(t)$ infected individuals, initially set to zero. This population evolves according to a continuous time branching process: in any given small time interval $dt$ $I_{k}(t)$ will increase by one with probability $\alpha_{k}I_{k}(t)dt$ and decrease by one with probability $\beta_{k}I_{k}(t)dt$. This represents infection and recovery ($\alpha$ and $\beta$ respectively). Parallel with this branching process, infected individuals may travel from one location, $i$, to a neighbouring location, $j$, with probability $\gamma T_{ij}I_{i}(t)dt$, resulting in an increment at location $j$ and a decrement at location $i$. We thus have what might be described as a “continuous time multitype branching process” [Kimmel and Axelrod, 2015].

We would like to determine the travel time from $a$ to $b$. Stated formally, suppose a novel diseases originates at time $t=0$, at node $a$. We would like to know the time at which the first infectious individual arrives in location $b$, that is to say the earliest time such that $I_{b}(t)>0$. We denote this epidemic arrival time (AT) as $\tau_{a}^{b}$. If no infected individuals ever arrives at $b$ we say that $\tau_{a}^{b}=\infty$. The distribution of $\tau_{a}^{b}$ will depend on the starting location $a$, the target location $b$, and the network parameters $\gamma,T$, $\alpha$ and $\beta$. We wish to understand this dependence and do so by considering the survival function ${S}_{a}^{b}(t)=\mathbb{P}(\tau_{a}^{b}>t)$.

Consider ‘Patient Zero,’ initially placed in location $a$ at time $t=0$. Now consider the state of the system a short time later, at time $t=dt$. With overwhelming probability, nothing will have happened during this small time window, and the survival times will still be governed by ${S}^{b}_{a}(t)$. With probability $dt\gamma T_{a,k}$ patient zero will have moved to location $k$, and the survival time distribution will instead be governed by ${S}^{b}_{k}(t)$. In the case where $k=b$, survival is now impossible, hence ${S}^{b}_{b}(t)=0$. Patient zero may also infect someone with probability $dt\alpha_{a}$. In this case, the probability of survival until time $t$ is given by the probability that neither traveler arrives at $b$. As both travelers are independent and identical, this probability is given by $[{S}^{b}_{a}(t)]^{2}$. Finally, it is possible that the initially infected patient zero may simply recover. This happens with probability $\beta_{a}$. In this case, no infected individual will ever reach $b$, and survival probability is $1$ for all time.

The probability of surviving for $t$ units of time after $dt$ is thus a superposition of the possible survival functions for the states the system could transition to. By definition the probability of surviving for $t$ units of time after time $dt$ is also precisely equal to ${S}^{b}_{a}(dt+t)$. Stated algebraically we thus have:

We refer to arrival times calculated using this method as ‘branching process arrival times’ (BP AT). These equations are equivalent to those given by Goldie & Coldman [Goldie and Coldman, 1979], who model the arrival of treatment resistant cancer cells via rare mutation (although in Goldie’s case, a much smaller collection of ‘types’ of individual are considered).

It is important to note that eq. (5) does not model the internal state of the system. While normal differential equations might be constructed by modelling the internal ‘state’ of the system at time $t$, and using this to model the system at $t+dt$, instead we merely observe that the survival function $S_{a}^{b}(dt+t)$ (whatever it might be), can be written as a superposition of the (as yet unknown) survival functions $S_{k}^{b}(t)$, and hence derive eq. (3). In so doing, we have lent heavily upon the assumption that system parameters are constant in time; that the behavior and survival curve associated with a single traveller observed at time $t$ is identical to the behavior and survival curve of a traveller observed at time $0$. For a discussion of how a similar approach might be taken in the context of time varying parameters, see Appendix A.

Equation 5, while accurate, is somewhat cumbersome and does not provide any closed form for values of interest such as expected arrival time and variance in arrival time. If we place no constraints on $\alpha$, $\beta$ and $\gamma T$ then it is easy to construct survival curves with complex topology such as fig. 4, and it can be shown that any valid survival curve (that is any curve of the form $0\leq f(t)\leq 1,f^{\prime}(t)\leq 0,f(0)=1$) can be well approximated, for suitably chosen network and parameter values (see appendix C). While mathematically interesting, such results are of little use to epidemiologists. More useful results can be obtained by placing some mild constraints on parameter values based on real world understanding.

The first assumption that we make is that we are not interested in epidemics that go extinct before spreading; we care only about those epidemics that reach international prevalence. For simplicity in the following discussion, this is achieved by assuming $\beta_{k}=0$ for all $k$. A discussion of non-zero $\beta_{k}$ can be found in Appendix D.

The second assumption that we make in what follows is that epidemic dynamics (infection and recovery) take place on on the time scale of days and weeks, while international travel is a rare occurrence, the vast majority of the human population taking international flights at most once or twice per year¹¹1With some individuals travelling significantly more often, and many others never flying at all.. Stated algebraically we have $\gamma T\ll\alpha$; travel is rarer than infection. This mimics the ‘rare mutation’ assumption made by Goldie & Coldman[Goldie and Coldman, 1979]when studying the the development of treatment resistant cancer cells, and rules out slow epidemics such as HIV.

We further assume that infection rate is equal across all locations: $\alpha_{k}=\alpha$. This greatly simplifies notation, and has limited impact upon our final conclusions.

Taken together, these three assumptions ($\beta_{k}=0,\alpha_{k}=\alpha,\gamma T_{ij}\ll\alpha$ are sufficient to avoid the more complex survival curve dynamics, as observed in fig. 4. Instead, we find $\dot{S}^{b}_{a}(t)\approx-\alpha{S}^{b}_{a}(t)\left(1-{S}^{b}_{a}(t)\right)$ and the arrival time distribution is well approximated by a logistic function:

Here $\mu_{a}^{b}$ is the mean arrival time (AT) for an epidemic starting in node $a$ and spreading to $b$. Determination of $\mu_{a}^{b}$ is of some considerable interest. This value is dependent upon the dynamics of the system when $(1-{S}^{b}_{a}(t))=O(\gamma)$; that is to say the region where the transport $\gamma TS^{b}_{a}(t)$ is no longer overwhelmed by logistic decay in probability.

So far we have explored the branching process model in what may be considered close to optimal conditions; we have assumed that $T$, $\alpha$ and $\gamma$ are exactly known, and that the susceptible population of each node is large enough so as not to limit epidemic growth and spread. Each of the above assumptions may be violated in one context or another, and for this reason it is critical to understand how far each can be stretched. Such understanding sheds light not only on the BP AT itself, but also on past distance metrics, which we have shown to rely on a similar theoretical foundations.

The first, and most likely assumption to be violated in practice is the notion that we know $\alpha$ and $\gamma$. While $\gamma$ can be determined to some reasonable level of accuracy via commercial flight data, $\alpha$ is a number that will vary from illness to illness and may be known only to a limited degree of accuracy, particularly in the early stages of a pandemic. Fortunately the model turns out to be robust against variations in both $\alpha$ and $\gamma$ – even when varying these parameters by a factor of five compared to the ‘true’ values used in simulations, predicted ATs remain highly correlated with those observed in simulation (see figure 7). When using incorrect $\alpha$ the constant of proportionality between $\tau_{i}$ and observed ATs is no longer equal to one, suggesting that the BP model robustly estimates $\alpha\tau_{i}$, such that the relative time of arrival at any two locations is largely independent of $\alpha$ and $\gamma$, but that incorrect estimates of $\alpha$ will lead to corresponding inaccuracy in the absolute values of $\tau_{i}$; if $\alpha$ is estimated a factor of two too high, then $\tau_{i}$ will be a factor of two too low.

Another systematic source of error is inaccuracies in flight network data. These can be produced due to uncounted or unregistered flights (or alternative forms of transportation), out of date data, or as a direct result of changes to individual flight plans as a result of the epidemic itself. In order to test robustness to noise in available flight data, we compare the results of full Gillespie simulation using flight network $T$, to ATs predicted using the perturbed matrix $\hat{F}$, where $\hat{T}_{ij}=T_{ij}e^{N(0,\xi)}$. Here $\xi$ is a constant representing noise level. As can be seen in Fig. 9, the resulting errors are modest.

While we have thus far considered arrival times in the context of the global aviation network, effective distance measures have also been considered on a much smaller scale in order to study the spread of antibiotic resistance between wards within a hospital [Friso Coerts, 2018]. In this case, assuming unbound growth becomes problematic, as the susceptible population in each node is rather small; a single ward might contain (for example) $10-30$ beds, a far cry from the tens or hundreds of thousands found when each node represents an entire city. When the local susceptible population is small, the epidemic is frequently no longer in an exponential growth phase by the time it spreads to a neighbouring node; ‘branching random walks’ are no longer independent because it is impossible for person $A$ to infect someone who has already been infected by person $B$. Population saturation inevitably leads to less accurate predictions, see figure 9.

In practice, for the SIR type models the parameter window where this concern is relevant is relatively narrow; infections which grow quickly enough to violate our assumptions are liable to burn through the entire local population and drive themselves to extinction before spreading between nodes. For SIS and SI type models, or any infection which is expected to reach a local endemic equilibrium, it is important to determine whether exponential growth is a good approximation before making use of the BP AT, or any distance metric that relies upon the same underlying unbound growth assumption.

Finally, while comparison to simulations provides an effective test bed, useful in terms of repeatability and certainty of data, we are, generally speaking more interested in the performance of models as they apply to the real world. Using flight data provided by Dirk Brockmann (private correspondence) we compare observed and predicted AT for both SARS and H1N1 outbreaks (see figure 10). We consider predictions made by both the branching process method described here (Eq. 5) as well as Brockmann & Helbing’s effective distance metric (Eq. 2, [Brockmann and Helbing, 2013]).

Ideally, in order to calculate the BP AT distribution, we would prefer to have access to the full transport matrix $T_{ij}$, the probability that an individual in location $i$ will travel to $j$ on any given day. Unfortunately such data is not available, and can not realistically be obtained. Instead, we have access to $F_{i,j}$; a matrix containing the total number of passengers travelling from $i$ to $j$ during a given time period. In order to approximate $T_{i,j}$ we normalize $F$, hence calculating the probability, $P_{i},j$, to fly from $i$ to $j$, conditioned on the assumption that we board some plane in airport $i$, and multiply by a constant $\gamma$, representing the probability of boarding some flight on a given day. In the real world, we might assume $\gamma$ to vary significantly from place to place, based on economic factors, prevalence of tourism, and the size of the population that any given airport is servicing. For the time being we approximate $\gamma$, and rely on the models low sensitivity to $\gamma$ to prevent difficulties.

We find that both BP AT and effective distance arrival time (‘ED AT’) methods give qualitatively similar results (see fig. 10), and that these results would appear in many cases plausible: for example, in the case of SARS, we observe that predicted arrival times in both South Korea and Hong Kong are low, while the arrival time in the USA is predicted to be marginally higher. This result is entirely plausible given airline traffic between the respective countries. In practice, predictions for both BP AT and ED AT do not match observed arrival times: SARS is reported to have arrived in the USA $54$ days after the first reports in China (11/16/2002), and in South Korea a full $160$ after these first reports. For SARS, Arrival times as reported from real world data correlate at best weakly with predictions(note the small $R^{2}<0.3$ in fig. 10 C,D), indicating either inaccuracy in the data or alternatively some complexity in the real world process not accounted for in current models. $R^{2}$ values are somewhat better when predicting H1N1 ($R^{2}\approx 0.55$).

We find no method (including ED) is able to reproduce the very high $R^{2}$ values reported in Brockmann & Helbing’s original paper, and that we are unable to recreate their figures (Figures 2D and 2E in the original paper), instead observing a far broader scatter. It is unclear if the observed discrepancies is due to differences in implementation of the algorithm, differences in the underlying data, or some other cause. While it is possible that inaccuracies in WAN data are to blame for these discrepancies, this explanation seems unlikely; as demonstrated in figure 9, arrival time predictions are generally stable, even to substantial changes in $F$. If we consider the effective distance metric (Eq 2), increasing $D_{i,j}$ by $5$, as is needed to correctly predict SARS late arrival time in South Korea, would require a $P_{i,j}$ two orders of magnitude smaller than our current value. Such drastic changes appear implausible.

Better understanding the spread of epidemics through the WAN allows for both real time forecasting (as has been used during the current COVID-19 pandemic), and also the possibility of network design, making changes to the WAN network so as to slow epidemic spread.

In studying the question of epidemic arrival time, a variety of models have been used, from the the most intricate agent based simulations[Broeck et al., 2011] to the intuitively appealing ‘distance’ based models [Brockmann and Helbing, 2013], and a number of analytic approaches in between [Chen et al., 2018, Iannelli et al., 2017, Gautreau et al., 2008]. While each of these models approaches the question of epidemic arrival time from a different view point, one common thread is the assumption of exponential growth; in the early stages of an epidemic, spread is governed by unbound epidemic growth, and rare transportation events. Population saturation, detailed viral dynamics, and a travellers tendency to return to their port of origin are ignored in all but the most detailed of models.

Given this basic premise, we were able to formulate the problem of epidemic arrival times in term of ‘Branching Processes’ [Kimmel and Axelrod, 2015, Goldie and Coldman, 1979], and calculate the full probability distribution of possible arrival times explicitly. These predictions match perfectly to corresponding simulations. If we further assume that air travel is rare and infection and recovery are common, it is possible to re-derive many of the results of past papers, and in some cases improve upon them. We are also able to give predictions in regions of parameter space where past methods are known to fail, and show that theoretical predictions provide a reasonable match to simulations, even when parameter values are altered significantly.

Unfortunately, when comparing to real world data we observe that the predictive powers of ED metrics may be significantly lower than previously reported [Brockmann and Helbing, 2013]. When compared to real world date, ED and BP like methods predict roughly 50% of variance in the cases of H1N1, and only 20% of variance when compared to SARS-2003, significantly lower than would be predicted given the modest intrinsic variance of our models. While the Branching Process Arrival Time makes very few assumptions and gives exact results given the underlying model assumed, these findings are suggestive of a gap in our knowledge- either in the data available for epidemic arrival times or (more likely) in the model itself. It seems likely that future efforts would be best directed towards determining what real world factors are currently unaccounted for, and which of these are most critical.

We gratefully acknowledge Dirk Brockmann for stimulating conversation, and for providing historical network and epidemic data. Full code for simulation based figures is available on github [Jamieson-Lane, 2020].