A Hybrid Model for Disease Spread and an Application to the SARS Pandemic

Pandemics can cause immense disruption and damage to communities and societies. Thus far, modeling of pandemics has focused on either large-scale difference equation models like the SIR and the SEIR models, or detailed micro-level simulations, which are harder to apply at a global scale. This paper introduces a hybrid model for pandemics considering both global and local spread of infections. We hypothesize that the spread of an infectious disease between regions is significantly influenced by global traffic patterns and the spread within a region is influenced by local conditions. Thus we model the spread of pandemics considering the connections between regions for the global spread of infection and population density based on the SEIR model for the local spread of infection. We validate our hybrid model by carrying out a simulation study for the spread of SARS pandemic of 2002-2003 using available data on population, population density, and traffic networks between different regions. While it is well-known that international relationships and global traffic patterns significantly influence the spread of pandemics, our results show that integrating these factors into relatively simple models can greatly improve the results of modeling disease spread.


1.1
Modeling of the spread of infectious disease typically falls into one of two categories. Analytically tractable models like the SEIR model are capable of capturing some globally important phenomena like the rate of spread of diseases using few parameters. However, they have a hard time reflecting differences in global spread due to local conditions. For example, it can be difficult to model different rates of spread in countries with different population densities and public health policies of variable strength and coordination. Network-or agent-based models are capable of reflecting details of individual conditions. However, modeling largescale global disease-spread using such models often runs into methodological problems like overfitting because of the vast number of possible parameters.

1.2
This paper proposes a granular, network-based hybrid model of disease spread in which individual regions are modeled as nodes in the network, and the spread of disease within nodes is modeled analytically (using a simplified derivative of the SEIR model) with the help of demographic parameters like population density. The properties of the network as a whole, like connectivity, are determined using real data on traffic between regions. We demonstrate the power of this approach by simulating the spread of SARS. One of the key takeaways is that the level of granularity has a significant effect on the success of networkor agent-based simulation models. For example, we show that modeling China as an individual node is unsuccessful, whereas breaking it up into constituent regions gives an impressive match to real infection data on SARS.

1.3
One of the great advantages of our model is its parsimony: it contains relatively few tweakable parameters compared with general agent-based models. At the same time it is capable of reproducing the important broad flows of disease. However, it is important to remember that exact reproduction of historical data is not the end-goal. Exceptions that do not correspond to real data provide insight into specific local phenomena that influence the progression of a pandemic, such as the actual timing of the first infected case in a country.

1.4
There is a vast literature on understanding the spread of disease using analytical and simulation models. In the next section we give a brief overview of the most common modeling methodologies, including differential equation models and simulation models, but here we discuss related research more generally. The most closely related to this work can be grouped into two categories. First, several researchers have simulated and analyzed the local spread of SARS in 2002Li et al 2004;Zhang et al 2005). In particular,  reproduce the situations in Singapore, Taipei, and Toronto individually, and compare with the actual transitions. This also ties in to a significant existing literature on local modeling of historical pandemics, like Influenza during the First World War (e.g. Chowell et al 2006;Massad et al 2006). Other examples also abound: Jenvald, et al (2007) use a virtual city based on Linköping, Sweden, considering the number of schools, age distribution, and household type; Longini, et al (2005) model population and contact processes based on Thailand census data, demographic information, and social network data; Kelso and colleagues model a real community in the southwest of Western Australia (Kelso et al 2009;Milne et al 2008).

1.5
The second category involves simulating global infection spread using international traffic data. For example, several papers use air travel data to estimate connectivity in a network (Ferguson et al 2006;Colizza et al 2007;Flahault et al 2006;Cooper et al 2006). However, these authors typically simulate a hypothetical global pandemic, with a focus on intervention policies; the focus of our research is to validate the simulation with real historical data.

1.6
Much existing research simulates infection in networks with reasonable properties, but not necessarily based on existing realworld data. For example, Bailey simulates epidemics in two dimensions, such as square grids (1965 (2002) also generate complex networks for simulation.

1.7
Another major theme of research has been on the effects of prevention and/or mitigation strategies. These typically compare a "base" simulation and an alternative simulation which considers some proposed strategy. For example,  use stochastic epidemic simulations to investigate the effectiveness of targeted antiviral prophylaxis to contain influenza. Kelso et al (2009) Figure 1 shows the spread of SARS as of April 8, 2003. In 20 of 29 countries/regions, 100% of total cases in the country were "imported" (as defined by WHO) from other countries (WHO Regional Office 2006).

1.9
The SARS pandemic is a particularly useful case study because we have high-fidelity data on the outbreak. First, the beginning and end of the pandemic are clear.

2.1
We first introduce the main existing methodologies used for modeling the spread of infectious disease before describing our approach in detail.

SIR Model
The classic SIR model, proposed by Kermack andMcKendrick in 1927 (Capasso 1993), posits three classes of agents; Susceptible, Infectious, and Removed. Susceptible agents (hereafter denoted S ) are vulnerable to a disease and have the potential to be infected. Infectious agents (I ) are currently infected and have the risk of infecting S . Removed (R ) agents are removed from the system -they are either dead or have acquired immunity.

2.2
Thus R-category agents are not infected again (R may also be called the Recovered category when we assume those agents are not dead). When R is not dead but has instead acquired immunity, the total population, (S + I + R ), is constant. The model assumes that agents in the set S are sometimes infected by a contact in I and change to R at a constant rate. This yields the expressions below for the transition of populations of these three classes. (1) where β is the rate of infection from S to I and λ is the rate of recovery from I to R . λ is inversely proportion to the average infectious period, τ : λ = τ -1 . When β/λ > 1 , the infection spreads since the susceptible become infected faster than the infected agents recover. The basic reproduction number, R 0 , is the average number of persons infected by a single infected person when the population has no immunity and no control against the infection (Yamamoto 2006). In the SIR differential equation model, the basic reproduction number is given by R 0 = Nβ/λ . If one infected person infects more than one susceptible person (i.e., R 0 > 1 ), secondary infection occurs and the infection spreads. On the other hand, if R 0 ≤1 , the disease converges.
Therefore R 0 = 1 is a threshold for spread.

SEIR Model
The SEIR model is a derivative of the SIR model. SIR does not consider the incubation period. Thus, when S is infected, it becomes I immediately and starts to infect other S (Wolfram Mathworld 2010). In the real world, there is some duration between the time that a person is infected and the time that he/she starts infecting others. The SEIR model denotes agents in the incubation period as belonging to class E (exposed) (Hethcote and Tudor 1980). The transition equation for the number of infected agents in the SIR model now works for the number of exposed agents, and the basic equations for the number of susceptible and recovered agents remain the same. A new equation has to be introduced for the number of infected agents that takes into account two different transitions: the transition from exposed to infected (increasing the number infected) and the transition from infected to recovered (decreasing the number infected). Now the transition equations become: (2)

2.3
Agent-based modeling provides an explicit, local method of understanding the spread of infection. It allows for fine-grained control over many aspects of the dynamic model of disease spread, including geographic factors and agent movements.  (2006a; 2006b) for simulating the spread of disease considering modules such as human activities, opportunity for contact between people in a society, disease state, and intervention to control the spread.

2.4
Network-based models typically represent agents as nodes on graphs and allow the connectivity structure of the graph to determine the possible spread of disease. For example, extending an SIR model to networks would involve allowing a susceptible vertex S to be infected by an infectious vertex I only if S is adjacent to I . Network-based models are useful in that they can reflect social and economic networks. People's behaviors and social contacts build the network and the infection route is on the network.

2.5
Our model uses local regions and interconnections between them. There are three possibilities for a new infection in a region; (1) infection from travelers from outside the region, (2) infection from returning travelers, and (3) infection from local persons. We denote infection types (1) and (2) as "global" infections and type (3) as "local" infections. Figure 2 shows the basic structure of our model. We note that while we are modeling stochastic spread of infection, the model itself outputs deterministic results which can be interpreted as a time series of the expected number of infections.

2.6
We assume that infection starts in a particular country or region and spreads from there. At each cycle, infections of all types can occur. Global infections (types 1 and 2) occur with frequencies that are dependent on the level of travel between regions, and local infections are mostly dependent on the population density of a region (details of the data used are below). Our local model is based on the concept behind the SEIR model. We consider the same four types of agents in each region: Susceptible, Exposed, Infectious, and Removed. When an infection occurs, agents are considered exposed. The model proceeds in time cycles t . The number of agents newly exposed in region i at time t through the global infection mechanism is modeled as EG i (t) = ∑ j I j (t) where the sum is over other regions, T ij is the sum of travelers from region i to region j and the number of travelers from region j to region i (since infection can occur through both arriving and returning travelers), and P G * (t) is a "global infection coefficient" at time t , described below.

2.7
Local infection follows a similar process, so that the number of agents newly exposed through the local infection mechanism at any time t is given by EL i (t) = S i (t)⋅I i (t)⋅P Li * (t) where P Li * (t) is a "local infection coefficient" (similar to the global infection coefficient, both are described in detail below).

2.8
It is assumed that agents go from exposed to infectious according to some incubation period that is disease-specific, and, similarly, from infectious to removed according to some disease-specific recovery period. For the purposes of this paper, we set these to 10 for both incubation period and infectious period, but these parameters can of course be varied for modeling other diseases.

2.9
As awareness of a disease spreads, it is likely that heightened awareness and prevention measures start to reduce the spread of infection. We model this in our global and local infection coefficients, by introducing a term that dampens the coefficient over where P G is a basic global infection coefficient, held constant across regions, and D G models the damping effect.

2.10
We use a similar equation for the local infection coefficient, P Li * (t) = P Li -(D L ⋅t) . In this case, D L is assumed constant across regions, but the basic local infection coefficient P Li is region-specific, and linear in the population density (assumed to be the main driver of local infection rates): P Li = ρ i ⋅C 1 + C 2 where ρ i is the population density of region i .

2.11
It is worth noting that the original SEIR model gives a similar type of equation for newly exposed agents E = β⋅S⋅I , where β is the infection rate. The main novelty here is the combination of modeling a declining infection rate, and treating each region separately.

3.1
There are several model parameters that need to be calibrated using real data. It is useful to consider some background information on the characteristics of SARS in this context.

3.2
The SARS Coronavirus causes general infection with Viremia, especially severe pneumonia and intestine infection. It is transmitted primarily through droplet infection. Due to its resistance to dryness, it can also be transmitted through air. It is thought that the incubation period of SARS is usually 2-10 days and the average is 5 days (Okada and Tashiro 2003). In the pandemic of 2002-2003, most countries reported a median incubation period of 4-5 days, and a mean of 4-6 days. In the incubation period, it is unlikely an infected person will spread the disease through droplet infection. The infectious period is thought to be about two weeks, with its peak from the 7th-10th day after infection (Okada and Tashiro 2003). Transmission efficiency appears to be greatest from severely ill patients who are experiencing rapid clinical deterioration, usually during the second week of illness. Maximum virus excretion from the respiratory tract occurs on about day 10 of illness and then declines to 0% by day 23. There are no reports on transmission beyond 10 days of fever resolution (University of Hong Kong 2010). The death rate varies by age group (SARS affects older patients much more severely), but the overall death rate was about 9.6% in the 2002-2003 SARS pandemic, significantly higher than that of seasonal Influenzas. Another notable feature of SARS is that it is believed that "superspreading" events, where a person infects many more than the average rate of infection, are a key component in its transmission. Our model does not deal explicitly with such levels of granularity, which may lead to some outlier predictions in areas where the law of large numbers does not take over. This is discussed further in Section 5.

3.3
It is thought that the origin of SARS was Guangdong in China, and the disease spread quickly from there to Hong Kong.

3.4
We hypothesize that population density of an area is positively correlated with the local infection rate, because higher population densities lead to more frequent contact. We test this hypothesis using data from Chinese provinces, Hong Kong, and Taiwan, the most significant infected regions. Figure 3 shows how the number of infections in different Chinese provinces varied greatly at the peak of the infection (from (Nations Online Project 2010)), which makes it necessary to treat the individual regions separately. Since 97% of infections occur in 6 provinces, we use data from these 6. They are Guangdong Province (the initial infected province), Beijing Municipality, Shanxi Province, Inner Mongolia Autonomous Region, Hebei Province, and Tianjin Municipality. Table 2 shows basic data on population and density for each of the provinces, Hong Kong, and Taiwan. Using these 6 provinces, Hong Kong, and Taiwan, we can reject the null hypothesis that there is no correlation between population density and infection rate at the 0.01 -level. Fitting C1 and C2 optimally is a difficult problem because of the complex dependence of infection rates on the multipliers. We manually tuned the parameters to provide a reasonable fit to the infection data from very early on in the infection cycle for different regions.    Table 4 shows the breakdown of total passenger traffic by type. We use land traffic data for the 6 provinces we are interested in, as shown in Table 5. We compute the share of each region in the total, where the total share is 100.

3.6
Then, based on Tables 4 and 5, we approximate the number of travelers between two regions by assuming that the share of a region is directly proportional to the number of travelers to the region. Also we assume that the share of passenger traffic by air is proportional to the share of passenger traffic by land.  We estimate travel between the different regions of China, Hong Kong, and Taiwan by using the share of the airport of each region in China in the national total. Table 6 shows the number of passengers using the main airport in 6 regions of China and the share in the national total. We apportion the number of travelers between China and Hong Kong or Taiwan according to the share. For example, the share of Beijing airport in the national total is 13.7859%.  For the preliminary experiment, we simulate with 6 regions in the Chinese mainland, Hong Kong, and Taiwan. The number of susceptible agents in each region/country is initially equal to the population of each country. Table 7 shows the summary for parameter values used in simulation. Note that these parameters are in terms of simulation cycles, which do not exactly correspond to real time. These simulation parameters were chosen to provide a good fit to data from this initial simulation, but we discuss below several inferences that can be made because many of the parameters are constant, exploiting the granularity of the model. Then, in the second part of this section, we use the same parameters to extend the model to 30 countries/regions, which provides a test for the parameters, allowing us to evaluate the benefits and drawbacks in a validation setting.  Table 2 T ij See Table 8 2.5 x 10 -7 Run_cycle 100 See Table 2 Incubation_period 10 Infection_period 10 P Li Depends on Density, C 1 , and C 2 5.0 x 10 -9 4.2 Figure 4 show the transition of the number of infected cases and the number of cumulative cases respectively, comparing real data and the results of our model. For the model we show data from time cycles 45 through 75. The results of China's 6 regions are summed up and the total is shown for China.
Density i D L Population i C 1 7.23 x 10 -9 C 2 7.69 x 10 -6 D G

4.3
The figure shows that our model captures both the dynamics of the spread of SARS, as well as the total numbers, very well. The peaks come in order: Hong Kong, China, and Taiwan. The model achieves this without using any special parameters that vary in different countries. Populations, densities, and travel data are all taken from the real world. The SARS epidemic started spreading from Hong Kong and immediately reached mainland China. The peak of Hong Kong comes earlier than that of China since the population density is higher. However, the curve decreases from some point because the percentage of susceptible agents in the population decreases and the percentage of Recovered agents increases. Then the number of Infectious agents decreases. After that, the number of infected agents in China increases. Because of its population, the number of infected agents at its peak in China is the largest among the three countries/regions. The peak in Taiwan is slightly delayed because of the time lag in the infection reaching Taiwan.

Region-wise Breakdown
4.4 Figure 5 shows the predicted (from the model) and actual number of cases for each of the eight modeled regions. While the fit is good for several of the most important (in terms of number of cases) regions, and therefore the overall numbers are good, there are some discrepancies for some of the regions that had relatively fewer cases. Specifically, the model underpredicts the number of cases for some of the less densely populated provinces of China (Shanxi, Inner Mongolia, and Hebei) and overpredicts for one of the more densely populated regions (Tianjin). There are idiosyncratic events associated with the spread of any pandemic, so it is not entirely surprising that some of the results do not match perfectly. The next section considers anomalies in more detail, where some data is available. But it is important to note that the level of granularity in modeling is very important. Figure 6 shows the difference in the model in two cases: one where the six infected provinces in China are modeled independently, and one where the six provinces are aggregated into one, using aggregated data on population density, travel etc. The figure clearly shows that the more granular model is a much better fit to the data.

4.5
As mentioned above, we use the parameters from the 8 region/country simulation to extend the model to 30 total countries (35 region/countries, since we continue to divide China into 6 regions). Again, we use real population, density, and international travel data from the 27 new countries (for Canada and Vietnam we use only Toronto and Hanoi, since only these regions had local spread cases (WHO Regional Office 2006)). Table 8 show the expected number of travelers between countries/regions. We again apportion the number of travelers between each region in mainland China and other countries based on the share of each region. 4.6 Figure 7 shows the infection route in our model. Most countries are infected from Hong Kong or Guangdong. Some countries are infected from other regions. For example, Vietnam, Mongolia, and Russia are infected from Beijing.  Figure 8 shows the comparison of the number of cumulative cases in simulation with real data. Especially for the significantly impacted countries, the number of cases corresponds well. In the real world, there were 8 countries/regions which had local spread. In the model, 18 countries/regions develop local spread. There are four true statistical outliers in the data in terms of number of cases predicted by the model versus number of cases experienced in reality. These are Singapore, Macao, Canada, and Japan.

Discussion of anomalies
We hypothesize that the outliers in this case are related to the nature of the spread of SARS. An early, chance outbreak, in a country or region can lead to significantly more cases than expected. Similarly, if a country manages to avoid a case of SARS for longer than predicted by international travel data, heightened awareness and prevention strategies will lead to many fewer cases than expected. For SARS this factor may be particularly important, because there is considerable evidence that some people infected with SARS are "super spreaders" who may affect the trajectory of spread. While an infected person infects, on average, 1-3 people (Okada and Tashiro 2003), some infected people pass the virus to many other people (WHO Regional Office 2006). Although it is not clear what causes someone to become a super spreader, it is suspected that a person who has a chronic illness such as diabetes is more likely to be a super spreader (Okada and Tashiro 2003

4.8
In each of the outlier cases, where the model makes a significantly different prediction than the actual trajectory of the pandemic, it turns our that the first reported case happened at a different time than would statistically be predicted by travel flows. While Macao, Japan, and Republic of Korea have large numbers of travelers from/to China and Hong Kong. these countries experienced much less infection than predicted by the model. It turns out that each of these countries experienced its first infection at a much later date than predicted, as shown in Table 9. The Republic of Korea first experienced an imported case on April 25th 2003 and Macao on May 5th 2003. Japan was never infected. These countries imported their first cases one or more April 25th 2003 and Macao on May 5th 2003. Japan was never infected. These countries imported their first cases one or more months after Vietnam, Canada, Taiwan, Singapore and the Philippines. Meanwhile, Canada, despite being less strongly linked by travel to China and Hong Kong, was infected on February 23rd, early in the pandemic (in fact, from the original super-spread event at the Metropole Hotel).

4.9
To provide some more weight to this hypothesis, we ran the model, but this time using the actual time of first infection in the country rather than travel data. Other than that, the parameters of the simulation remained the same. Figure 9 shows that the cumulative number of cases from the model then correspond better to real data. Here we discuss two of them, and how they might affect the results. First, if we look at data from seasonal flu cases, we find that Canada typically has a large number of cases, and the United States has the largest number of influenza isolates (World Health Organization 2010). Both of these suggest that the local infection coefficient may be higher in Canada and the United States than other countries. Indeed, this could have been an additional factor in the surprisingly large number of Canadian cases. However, the United States was surprising, because, although it imported 27 cases, the infection did not spread locally. This may indicate that the quarantining and isolation measures employed worked effectively.

4.11
A second interesting point is that Singapore and Vietnam both report many fewer cases than predicted by the model. This may be partly explained by their lower propensity to spread infection, again as evidenced by seasonal flu data. There may also have been a significantly bigger push to hospitalize and keep patients confined, weakly evidenced by the fact that the proportion of those infected who were healthcare workers in these two countries (41% and 57% in Singapore and Vietnam respectively) was much higher than other countries (21%).

5.1
We have discussed a hybrid network and local model for the spread of pandemics, and applied it to the case of SARS. When parameters are calibrated to real data on populations, densities, and traffic, we show that the model reproduces many of the key dynamics of the spread of SARS in 2002 and 2003, while remaining parsimonious, and therefore useful for understanding the root causes of why pandemics spread in the way they do. Both the successes and the failures of the simple model provide insights into pandemic spread. For example, it is clear that it is important to model international traffic to understand the pathways of spread. At the same time, for any particular pandemic, individual idiosyncrasies can come into play. For example, the importance of super-spreaders in SARS is reflected in the fact that the time of first infection in a country plays a big role in how many people get infected. The other major takeaway from this work is that the level of granularity in the network structure of the model has a significant impact. For example, treating China as one large entity leads to poorer prediction, but at the same time specializing all the way down to cities would end up requiring too much data to accurately calibrate the model, without significantly improving the predictivity of the model.

5.2
Our model is quite general. It should be applicable without major changes to many different countries or regions, and even to different granularities of population centers, like cities. It should also apply to epidemics with different properties than SARS. We also believe that the idea of using hybrid network-based simulation models has great potential in interventional settings, because it is easy to experiment with how changing certain parameters like traffic flow between different parts of the network, affects the global outcomes for different types of epidemic phenomena.

Notes
A Java framework for our simulation, along with instructions for how to use it, can be found at http://www.cs.rpi.edu/~moorthy/SARS-Program/.