FLUed: A Novel Four-Layer Model for Simulating Epidemic Dynamics and Assessing Intervention Policies

Abstract

From the 2003 severe acute respiratory syndrome (SARS) epidemic, to the 2009 swine-origin influenza A (H1N1) pandemic, to the projected highly pathogenic avian influenza A event, emerging infectious diseases highlight the importance of computational epidemiology to assess potential intervention policies. Hence, an important and timely research goal is a general-purpose and extendable simulation model that integrates two major epidemiological factors—age group and population movement—and substantial amounts of demographic, geographic, and epidemiologic data. In this paper, we describe a model that we have named FLUed for Four-layer Universal Epidemic Dynamics that integrates complex daily commuting network data into multiple age-structured compartmental models. FLUed has four contact structures for simulating the epidemic dynamics of emerging infectious diseases, assessing the potential efficacies of various intervention policies, and identifying the potential impacts of spatial-temporal epidemic trends on specific populations. We used data from the seasonal influenza A and 2009 swine-origin influenza A (H1N1) epidemics to validate model reliability and suitability and to assess the potential impacts of intervention policies and variation in initial outbreak areas for novel/seasonal influenza A in Taiwan. We believe that the FLUed model represents an effective tool for public health agencies responsible for initiating early responses to potential pandemics.

1. Introduction

After emerging in Mexico in April 2009, the swine-origin influenza A (H1N1) virus rapidly spread to more than 200 countries, causing over 18,449 laboratory-confirmed deaths August 1, 2010 [1]. According to the World Health Organization (WHO) [2], up to two billion people may be susceptible to the next highly pathogenic avian influenza A (H5N1) virus, and the mortality rate could be as high as 65%. A novel influenza virus is inevitable yet impossible to predict, and health officials can only guess where, when, and in what form the virus will emerge, or how it might threaten the health of individuals [3]. From the 2003 severe acute respiratory syndrome (SARS) epidemic, to the 2009 swine-origin influenza A (H1N1) pandemic, to the projected highly pathogenic avian influenza A (H5N1, H7N9) event, emerging infectious diseases highlight the importance of computational epidemiology for assessing potential intervention policies [4, 5]. Thus, epidemic simulation models that integrate epidemiological factors with substantial amounts of statistical data represent a timely research topic attracting considerable resources.

To date, public health officials and computational epidemiologists have focused on two modeling factors: age group and population movement [6–11]. Since individuals in the same-age group tend to have similar postinfection symptoms and epidemic characteristics, many researchers use this factor to capture heterogeneity in epidemic models of emerging infectious diseases, to assess the potential efficacies of various intervention policies, and to identify the potential impacts of temporal epidemic trends on specific populations [6–12]. For example, the Mexican population segment that was most affected by the 2009 swine-origin influenza A (H1N1) virus consisted of youth below the age of 15—of all individuals affected by the first infection wave, 61% were children and 29% adults [13, 14]. By contrast, SARS symptoms in children during the 2003 outbreak were mild and short-lived [15]. Results from serological analyses indicate that levels of antibodies against the SARS coronavirus were higher in children than those in adults in both infected and healthy children [16].

Population movement is another important factor in modeling epidemic dynamics and spreading situations. Over the past three decades, many countries have experienced a rapid increase in the number of commuters for work and other purposes, especially among young adults [17]—a phenomenon perceived as supporting the spread of viruses over long distances [18]. Commuting is marked by strong spatial-temporal regularity; regardless of travel distance or time, most commuters stick to repetitive patterns that allow researchers to study the epidemic dynamics and clustering effects of emerging infectious diseases [5, 19, 20], as well as intervention policy activation timing and targeting [21–24]. For example, Viboud et al. [25] analyzed seasonal influenza data for various American states from 1972 to 2002 to determine the influence of interstate commuting on seasonal influenza epidemic dynamics and found a stronger correlation between regional infection spread and rate of residence-workplace movement than that of geographic distance. In other words, higher commuting frequency increases the likelihood of a larger number of seasonal influenza cases for any given year. According to their results, if the initial case occurs in a state that has a large population and large volume of outbound traffic, then the seasonal influenza will quickly affect all the fifty states, but if it occurs in a sparsely populated state with a small volume of outbound traffic, the potential for spreading to all states will be significantly less. In addition to providing a foundation for empirically verifying the ways that commuting networks affect seasonal influenza epidemic dynamics, Viboud et al.’s findings also verify that transportation networks and commuting populations strongly influence seasonal influenza outbreak scale and expansion speed.

The two most commonly used methods for modeling epidemic dynamics are network based and compartmental based. Each has its own limitations [22, 26, 27]. A compartmental model is suitable for discussing dynamic variation across individuals in the same compartment but weak in terms of modeling individual heterogeneity and addressing human commuting patterns. Any two individuals in a compartment are assumed to have direct contact and interaction, which is not true in the real world. Furthermore, since movement and activity are location dependent, phenomena cannot be simulated by a compartmental model that assumes a homogeneous population distribution. By contrast, network models are appropriate for introducing individual heterogeneity, but they are computation intensive and time-consuming when simulating the behaviors of individuals with multiple attributes in large-scale social environments. Many efforts have been made to match individual and population behaviors with heterogeneity and computation requirements when simulating epidemic dynamics and assessing intervention policies [28–30].

This paper has four objectives: (a) to describe a general-purpose¹ and extendable² model that integrates multiple age-structured compartmental models into complex daily commuting network data for simulating the epidemic dynamics and spreading situations of newly emerging infectious diseases; (b) to use data from the 2008-2009 seasonal influenza A and 2009 swine-origin influenza A (H1N1) epidemics to estimate model parameters; (c) to assess the potential impacts of different intervention policies and initial outbreak areas on the Taiwan-wide dynamics of the 2009 swine-origin influenza A (H1N1) epidemic, including intervention activation time and decreased transmission rate; and (d) to apply two quantitative indicators to assess the performances of different intervention policies. Our results not only support the use of the proposed four-layer framework for computational epidemiology applications but also support the framework’s potential utility for intervention policy assessments.

2. The Model

To perform detailed analyses of the specific epidemic dynamics and spreading situations of emerging infectious diseases in individual countries, we developed a general-purpose and extendable model that we named Four-layer Universal Epidemic Dynamics, or FLUed.³ FLUed uses a mix of multiple age-structured compartmental models and daily commuting networks to capture complex demographic, geographic, and biological properties, including population movement and epidemiological progression. As shown in Figure 1, layer 1 individuals in a defined location are organized according to age; a compartmental model is used to simulate the epidemic dynamics of each age group. A location can be defined as a region, town, city, county, or state, but all the locations must be of the same type. The layer 2 focus is on contact patterns among age groups in the same location. Population density and commuting volume between any two locations are added to layer 3 to study the effects of crossregional interaction. A geographic information system (GIS) user interface is integrated into layer 4 to visualize daily national commuting networks, with nodes⁴ representing locations and links representing daily commutes between them.

Figure 1 

FLUed model framework. Infection information is highest in layer 1 and lowest in layer 4 (opposite of location information).

The model’s concept framework consists of simulation flow (blue box), transportation, and census databases respectively, obtained from the Republic of China (ROC) Transportation Institute and Interior Ministry (red box) and data flow relationships (green box) among the four FLUed model layers, internal data structures, and external databases (Figure 2) [34]. Most data sets consist of spatial locations and census information, which are used to establish geographic and demographic categories; each model parameter belongs to at least one of the two. Since our layer 1 focus is on epidemiological progression at an individual level, we used standard expert-based parameters⁵ in compartmental models associated with emerging infectious diseases instead of the transportation or census databases [21]. In layer 2, the numbers of individuals in each location and percentages of individuals in each age group were determined from census data. Layer 3 information on the numbers of individuals commuting between locations on a daily basis were taken from the R.O.C. transportation database. The layer 4 national commuting network was constructed with information from the same transportation database and a GIS.

Figure 2 

FLUed model implementation framework.

Two methods were used to input and alter model parameters: a GUI simulation console (Figure 3) and input files (i.e., Excel spreadsheet files containing data for the scenario being analyzed). Default values for each model parameter were automatically loaded at startup. Model parameters initialized at the beginning of each simulation were (a) initial outbreak conditions, including the name of the location and number of infected persons in an age group identified by the surveillance system; (b) epidemiological parameters at different layers, including transmission, latent, and removal rates according to an age-structured compartmental model for each age group, contact rates between age groups, and regional contact probabilities between locations; and (c) GIS output maps and time-series charts for locations of interest and severity indicators to be monitored (e.g., daily infected cases, daily new cases, and epidemic velocity⁶ and acceleration rate⁷). Daily epidemiological progressions were monitored in terms of the sizes and locations of red dots on a GIS map, curves on time-series charts, and output panels showing numbers of infected individuals at different times in different locations.

Figure 3 

FLUed model simulation console.

We used a modified Susceptible-Latent-Infectious-Removed (SLIR) compartmental model to represent different infection stages among individuals in the same-age group in the same location. Individual epidemiological status was initially set at Susceptible (vulnerable to infection but not yet infected), followed by Latent (infected but unable to infect others), Infectious (capable of infecting other individuals), and Removed (i.e., recovered, deceased, or otherwise not posing any further threat). The numbers of pathogens carried by Susceptible-to-Latent hosts were initially insufficient for active transmission to other Susceptible hosts but eventually reached levels causing hosts to become Infectious to begin infecting other Susceptible hosts and to move toward a Removed status.

The dynamics of the four epidemiological statuses over time are expressed as (1a)–(1d). They have five features: (a) discrete time interval 𝑑𝑡 is assumed to be one day; (b) at time 𝑡, the population of interest is divided into four compartments—𝑆(𝑡), 𝐿(𝑡), 𝐼(𝑡), and 𝑅(𝑡)—corresponding to the four epidemiological statuses; since our SLIR model is a closed system, 𝑆(𝑡) + 𝐿(𝑡) + 𝐼(𝑡) + 𝑅(𝑡) = 𝑁, with 𝑁 a constant representing the entire population; (c) transmission rate 𝛽 is a constant representing how fast Susceptible individuals become infected and acquire a Latent status; (d) latent rate 𝜃 is a constant used to determine transformation speed from Latent to Infected; and (e) removal rate 𝛼 is a constant used to determine transformation speed from Infected to Recovered. Consider𝑑𝑆(𝑡)(𝑡)𝐼(𝑡)𝑑𝑡=−𝛽𝑆𝑁,𝑑𝐿(𝑡)(𝑡)+(𝑡)𝐼(𝑡)𝑑𝑡=−𝜃𝐿𝛽𝑆𝑁,𝑑𝐼(𝑡)(𝑡)(𝑡),(𝑡)𝑑𝑡=−𝛼𝐼+𝜃𝐿𝑑𝑅(𝑡).𝑑𝑡=𝛼𝐼 (1a)(1b)(1c)(1d)

A system dynamics flowchart for the first layer of the FLUed model is shown in Figure 4. Note our modification in consideration of self-motivated hospitalization (i.e., individuals who visit hospitals or clinics during infectious disease outbreaks regardless of their infection status). Depending on diagnostic accuracy, some self-motivated individuals are confirmed as Infectious and receive medical treatment in advance, thus altering transmission and removal rates for certain populations. We added three features to integrate this factor into the model: (a) a real number investigation constant 𝑠 (0≤𝑠≤1.0) representing the percentage of a population that goes to a hospital or clinic before becoming ill; (b) a real number detection constant 𝑐 (0≤𝑐≤1.0) used to determine the percentage of a population confirmed as Infectious; and (c) a time delay constant 𝑇, indicating the amount of time between a patient with symptoms visiting a hospital or clinic and the time an infection is confirmed. The default values of parameters 𝑠 and 𝑐 are both 0.6, meaning that 60% of the infected population is prone to visiting hospitals or clinics for medical advice, and 60% of those visitors are correctly diagnosed as carrying the pathogen. The default value of parameter 𝑇 is 3, meaning that it takes three days to confirm that a hospital or clinic patient with symptoms is infected. In our simulations, the numbers of correctly diagnosed patients were equal to the numbers of confirmed real-world cases identified by emerging infectious disease surveillance systems.

Figure 4 

System dynamics flowchart: modified SLIR compartmental model in layer 1 of the FLUed model.

In consideration of preventive health care efforts, we added a feature in which 𝐿 status individuals are moved to either an 𝐼₁ (infected and prone to visiting hospitals or clinics for medical advice) or 𝐼₂ (infected but not prone) status based on whether or not they visit a hospital or clinic; this feature is expressed as investigation proportion 𝑠. 𝐼₁ individuals are identified as either 𝐼₁₁ (correctly diagnosed as carrying the pathogen) or 𝐼₁₂ (incorrectly diagnosed—in other words, false negatives); this feature is expressed as detection proportion 𝑐. Note that regardless of positive or negative diagnoses, a 𝑇 period of time must elapse prior to confirmation. For individuals who are correctly diagnosed with influenza, some medical treatments and public health practices can control symptoms, prevent further complications, and stop them from spreading the virus to their families, friends, classmates, and coworkers. Correctly diagnosed patients in simulations are the equivalent of confirmed cases in real-world influenza surveillance systems. Thus, the difference between 𝐼₁₁ and either 𝐼₂, or 𝐼₁₂ statuses is the transmission rate. 𝐼₁₁, 𝐼₁₂ and 𝐼₂ all eventually change to status 𝑅. The extended SLIR model can be expressed as (2a)–(2h). Consider𝑑𝑆(𝑡)𝑆(𝑡)𝛽𝑑𝑡=−²𝐼₂(𝑡)+𝛽₁₁𝐼₁₁(𝑡)+𝛽₁₂𝐼₁₂(𝑡)𝑁,𝑑𝐿(𝑡)(𝑡)+𝑆(𝑡)𝛽𝑑𝑡=−𝜃𝐿²𝐼₂(𝑡)+𝛽₁₁𝐼₁₁(𝑡)+𝛽₁₂𝐼₁₂(𝑡)𝑁,𝑑𝐼₀(𝑡)𝑑𝑡=−𝐼₀(𝑡)+𝜃𝐿(𝑡),𝑑𝐼₁(𝑡)𝑑𝑡=−𝐼₁(𝑡−𝑇)+𝑠𝐼₀(𝑡),𝑑𝐼₂(𝑡)𝑑𝑡=−𝛼₂𝐼₂(𝑡)+(1−𝑠)𝐼₀(𝑡),𝑑𝐼₁₁(𝑡)𝑑𝑡=−𝛼₁₁𝐼₁₁(𝑡)+𝑐𝐼₁(𝑡−𝑇),𝑑𝐼₁₂(𝑡)𝑑𝑡=−𝛼₁₂𝐼₁₂(𝑡)+(1−𝑐)𝐼₁(),(𝑡)𝑡−𝑇𝑑𝑅𝑑𝑡=𝛼₁₁𝐼₁₁(𝑡)+𝛼₁₂𝐼₁₂(𝑡)+𝛼₂𝐼₂(𝑡). (2a)(2b)(2c)(2d)(2e)(2f)(2g)(2h)

Depending on age range, individual infection properties differ in terms of epidemiological parameters such as transmission and removal rates. We considered two age-related features: the transmission rates 𝛽_11𝑝𝑞, 𝛽_12𝑝𝑞, and 𝛽_2𝑝𝑞, which represent cross-age group infections and the relative age group percentage 𝜒_𝑝, which affects the potential for cross-age infections. To distinguish among parameters for individuals in different age groups, we also added a subscript to each (2a)–(2h) parameter (with the exception of 𝑇)—for example, we changed parameter 𝑆(𝑡) to 𝑆_𝑝(𝑡) for age group 𝑝. We assumed three age groups: children (birth to 14), adults (15 to 64), and seniors (65 and older). Transmission rates between age groups were differentiated to capture the complexity of infections across them. We added two subscripts to transmission rate 𝛽 to create 𝛽_𝑝𝑞: 𝑝 for the age of an infectious individual, and 𝑞 for the age of the individual being infected (Figure 5). Epidemiological parameters used in the first layer were also used to model infections across age groups.

Figure 5 

System dynamics flowchart between two modified age-structured SLIR compartmental models in layer 2 of the FLUed model. Red lines indicate parameters for other age groups. Blue curves indicate relative percentages of each age group within the total population at each location.

To construct layer 2 of the FLUed model, we revised (2a) and (2b) to (3a) and (3b), respectively, without making any other changes to the (2a)–(2h) subequations. Consider𝑑𝑆_𝑝𝑆𝑑𝑡=−^𝑝𝜒_𝑝𝜒_𝑝𝛽_2𝑝𝑝𝐼_2𝑝+𝛽_11𝑝𝑝𝐼_11𝑝+𝛽_12𝑝𝑝𝐼_12𝑝𝑁_𝑝−𝑆_𝑝∑_𝑞≠𝑝𝜒_𝑞𝜒_𝑝𝛽_2𝑞𝑝𝐼_2𝑞+𝛽_11𝑞𝑝𝐼_11𝑞+𝛽_12𝑞𝑝𝐼_12𝑞𝑁_𝑞,𝑑𝐿_𝑝𝜃𝑑𝑡=−^𝑝𝐿_𝑝+𝑆_𝑝𝜒_𝑝𝜒_𝑝𝛽_2𝑝𝑝𝐼_2𝑝+𝛽_11𝑝𝑝𝐼_11𝑝+𝛽_12𝑝𝑝𝐼_12𝑝𝑁_𝑝+𝑆_𝑝∑_𝑞≠𝑝𝜒_𝑞𝜒_𝑝𝛽_2𝑞𝑝𝐼_2𝑞+𝛽_11𝑞𝑝𝐼_11𝑞+𝛽_12𝑞𝑝𝐼_12𝑞𝑁_𝑞. (3a)(3b)

In layer 3 we focused on the impacts of daily commuting between locations on the epidemic dynamics of emerging infectious diseases. This layer reflects three assumptions⁸: (a) emerging infectious viruses are transmitted via airborne droplets, and commuters are capable of infecting other individuals along their standard routes; (b) according to census statistics obtained from the ROC Interior Ministry, hosts commute over longer distances than individuals who stay at home or travel to local centers such as schools and tend to come into contact with individuals in the same-age group along their routes and at their destinations; and (c) higher contact frequencies exist among individuals in more densely populated areas. Accordingly, layer 3 consists of four features associated with daily commuting: (a) 𝜎(𝑝), a binary value denoting whether age group 𝑝 is the commuter age group (our assumption was that children and seniors are less likely than adults to commute on a daily basis, making adults the most likely carriers of pathogens between locations); (b) 𝑤_𝑗,𝑖, indicating how many individuals commute from location 𝑗 to location 𝑖 on a daily basis; (c) 𝜂^𝑖, a weighting factor representing the average number of contacts among individuals in location 𝑖 on a daily basis; and (d) 𝑑^𝑖, a relative density representing the location 𝑖 population as a percentage of the population of the largest location in the commuting network.

For all 𝑖 and 𝑗 locations in a commuting network 𝑤, we used a geodemographic weight equation, 𝜎(𝑝)𝑑^𝑖𝜂^𝑖(𝑤_𝑗,𝑖/∑_𝑘≠𝑗𝑤_𝑗,𝑘), to measure the effects of commuting on 𝑖 and 𝑗 population interactions. The term 𝑤_𝑗,𝑖/∑_𝑘≠𝑗𝑤_𝑗𝑘 is the ratio of commuters between locations 𝑗 and 𝑖 to commuters between 𝑗 and all other locations. If location 𝑖 is a large urban center, 𝑤_𝑗,𝑖/∑_𝑘≠𝑗𝑤_𝑗,𝑘 will be large; if 𝑖 is a suburb or rural location, it will be small. Public health policies involving transportation can be tested by changing contact rates among population centers in the third layer.

Figure 6 presents the system dynamics flowchart for the FLUed model’s third layer; to construct it we, respectively, revised (3a) and (3b) to (4a) and (4b). Note the addition of a geodemographic weight factor on the third line of each equation. All other (2a)–(2h) subequations are identical. Consider𝑑𝑆^𝑖_𝑝=𝑑𝑡−𝑆^𝑖_𝑝𝑑^𝑖𝜒^𝑖_𝑝𝜒^𝑖_𝑝𝛽^𝑖_2𝑝𝑝𝐼^𝑖_2𝑝+𝛽^𝑖_11𝑝𝑝𝐼^𝑖_11𝑝+𝛽^𝑖_12𝑝𝑝𝐼^𝑖_12𝑝𝑁^𝑖_𝑝⁻¹−𝑆^𝑖_{𝑝𝑞≠𝑝}𝑑^𝑖𝜒^𝑖_𝑞𝜒^𝑖_𝑝𝛽^𝑖_2𝑞𝑝𝐼^𝑖_2𝑞+𝛽^𝑖_11𝑞𝑝𝐼^𝑖_11𝑞+𝛽^𝑖_12𝑞𝑝𝐼^𝑖_12𝑞𝑁^𝑖_𝑞⁻¹−_{(𝑖)𝑗∈𝑁}𝜎𝑝𝑑^𝑖𝜂^𝑖𝑆^𝑖_𝑝𝑤_𝑗,𝑖∑_𝑘≠𝑗𝑤_𝑗,𝑘×𝛽^𝑗_2𝑝𝑝𝐼^𝑗_2𝑝+𝛽^𝑗_11𝑝𝑝𝐼^𝑗_11𝑝+𝛽^𝑗_12𝑝𝑝𝐼^𝑗_12𝑝𝑁^𝑗_𝑝⁻¹,𝑑𝐿^𝑖_𝑝=𝑑𝑡−𝜃^𝑖_𝑝𝐿^𝑖_𝑝+𝑆^𝑖_𝑝𝑑^𝑖𝜒^𝑖_𝑝𝜒^𝑖_𝑝𝛽^𝑖_2𝑝𝑝𝐼^𝑖_2𝑝+𝛽^𝑖_11𝑝𝑝𝐼^𝑖_11𝑝+𝛽^𝑖_12𝑝𝑝𝐼^𝑖_12𝑝𝑁^𝑖_𝑝⁻¹+𝑆^𝑖_{𝑝𝑞≠𝑝}𝑑^𝑖𝜒^𝑖_𝑞𝜒^𝑖_𝑝𝛽^𝑖_2𝑞𝑝𝐼^𝑖_2𝑞+𝛽^𝑖_11𝑞𝑝𝐼^𝑖_11𝑞+𝛽^𝑖_12𝑞𝑝𝐼^𝑖_12𝑞𝑁^𝑖_𝑞⁻¹+_{𝑗∈𝑁(𝑖)}𝜎𝑝𝑑^𝑖𝜂^𝑖𝑆^𝑖_𝑝𝑤_𝑗,𝑖∑_𝑘≠𝑗𝑤_𝑗,𝑘×𝛽^𝑗_2𝑝𝑝𝐼^𝑗_2𝑝+𝛽^𝑗_11𝑝𝑝𝐼^𝑗_11𝑝+𝛽^𝑗_12𝑝𝑝𝐼^𝑗_12𝑝𝑁^𝑗_𝑝⁻¹. (4a)(4b)

Figure 6 

System dynamics flowchart of modified age-structured SLIR compartmental model in layer 2 of the FLUed model. Properties associated with commuting between two locations are indicated by red lines. Additional location properties are indicated by blue lines.

Based on national short- and long-distance commuting data, network nodes can represent individual towns and link attributes can represent geodemographic weights and commuting volumes between any two towns. As shown in Figure 7, layer 4 of the FLUed model consists of 409 towns and 19,014 links that can be manipulated to test the effects of various transport intervention policies, population movement strategies, and commuting restrictions. After combining the four layers, the model can be expressed as (5a)–(5h). Model parameters and their default values are listed in Table 1. The simulation time complexity of the model is 𝑂(𝛼²𝑛𝑘𝜏), a low-order polynomial function in which 𝛼 denotes the number of age groups, 𝑛 the number of commuting network nodes, 𝑘 the average degree of the daily commuting network, and 𝜏 the total number of time steps in an epidemic simulation experiment. One has the following:𝑑𝑆^𝑖_𝑝=𝑑𝑡−𝑆^𝑖_𝑝𝑑^𝑖𝜒^𝑖_𝑝𝜒^𝑖_𝑝𝛽^𝑖_2𝑝𝑝𝐼^𝑖_2𝑝+𝛽^𝑖_11𝑝𝑝𝐼^𝑖_11𝑝+𝛽^𝑖_12𝑝𝑝𝐼^𝑖_12𝑝𝑁^𝑖_𝑝⁻¹−𝑆^𝑖_{𝑝𝑞≠𝑝}𝑑^𝑖𝜒^𝑖_𝑞𝜒^𝑖_𝑝𝛽^𝑖_2𝑞𝑝𝐼^𝑖_2𝑞+𝛽^𝑖_11𝑞𝑝𝐼^𝑖_11𝑞+𝛽^𝑖_12𝑞𝑝𝐼^𝑖_12𝑞𝑁^𝑖_𝑞⁻¹−𝜎𝑝𝑑^𝑖𝜂^𝑖𝑆^𝑖_{𝑝𝑗∈𝑁(𝑖)}𝑤_𝑗,𝑖∑_𝑘≠𝑗𝑤_𝑗,𝑘×𝛽^𝑗_2𝑝𝑝𝐼^𝑗_2𝑝+𝛽^𝑗_11𝑝𝑝𝐼^𝑗_11𝑝+𝛽^𝑗_12𝑝𝑝𝐼^𝑗_12𝑝𝑁^𝑗_𝑝⁻¹,𝑑𝐿^𝑖_𝑝𝜃𝑑𝑡=−^𝑖_𝑝𝐿^𝑖_𝑝+𝑆^𝑖_𝑝𝑑^𝑖𝜒^𝑖_𝑝𝜒^𝑖_𝑝𝛽^𝑖_2𝑝𝑝𝐼^𝑖_2𝑝+𝛽^𝑖_11𝑝𝑝𝐼^𝑖_11𝑝+𝛽^𝑖_12𝑝𝑝𝐼^𝑖_12𝑝𝑁^𝑖_𝑝⁻¹+𝑆^𝑖_{𝑝𝑞≠𝑝}𝑑^𝑖𝜒^𝑖_𝑞𝜒^𝑖_𝑝𝛽^𝑖_2𝑞𝑝𝐼^𝑖_2𝑞+𝛽^𝑖_11𝑞𝑝𝐼^𝑖_11𝑞+𝛽^𝑖_12𝑞𝑝𝐼^𝑖_12𝑞𝑁^𝑖_𝑞⁻¹+𝜎𝑝𝑑^𝑖𝜂^𝑖𝑆^𝑖_{𝑝𝑗∈𝑁(𝑖)}𝑤_𝑗,𝑖∑_𝑘≠𝑗𝑤_𝑗,𝑘×𝛽^𝑗_2𝑝𝑝𝐼^𝑗_2𝑝+𝛽^𝑗_11𝑝𝑝𝐼^𝑗_11𝑝+𝛽^𝑗_12𝑝𝑝𝐼^𝑗_12𝑝𝑁^𝑗_𝑝⁻¹,𝑑𝐼^𝑖_0𝑝𝑑𝑡=−𝐼^𝑖_0𝑝+𝜃^𝑖_𝑝𝐿^𝑖_𝑝,𝑑𝐼^𝑖_1𝑝𝑑𝑡=−𝐼^𝑖_1𝑝(𝑡−𝑇)+𝑠^𝑖_𝑝𝐼^𝑖_0𝑝,𝑑𝐼^𝑖_2𝑝𝑑𝑡=−𝛼^𝑖_2𝑝𝐼^𝑖_2𝑝+1−𝑠^𝑖_𝑝𝐼^𝑖_0𝑝,𝑑𝐼^𝑖_11𝑝𝑑𝑡=−𝛼^𝑖_11𝑝𝐼^𝑖_11𝑝+𝑐^𝑖_𝑝𝐼^𝑖_1𝑝(𝑡−𝑇),𝑑𝐼^𝑖_12𝑝𝑑𝑡=−𝛼^𝑖_12𝑝𝐼^𝑖_12𝑝+1−𝑐^𝑖_𝑝𝐼^𝑖_1𝑝(𝑡−𝑇),𝑑𝑅^𝑖_𝑝𝑑𝑡=𝛼^𝑖_11𝑝𝐼^𝑖_11𝑝+𝛼^𝑖_12𝑝𝐼^𝑖_12𝑝+𝛼^𝑖_2𝑝𝐼^𝑖_2𝑝. (5a)(5b)(5c)(5d)(5e)(5f)(5g)(5h)

Table 1 

FLUed model expert-identified epidemic parameters used for fitting simulation curves with actual 2008-2009 seasonal influenza A and 2009 swine-origin influenza A (H1N1) curves in Taiwan.

Figure 7 

National commuting network in Taiwan.

3. Model Validation and Intervention Policy Evaluation

We used case data confirmed by the Taiwan Centers for Disease Control (TCDC) for 2008-2009 seasonal influenza A and 2009 swine-origin influenza A (H1N1) to estimate FLUed model parameters and to test model reliability and generalizability. We calibrated parameters for both influenzas to create small ranges based on parameters normally used in extended compartmental model settings [16]. Parameter and value summaries are given in Table 1. The transmission rates 𝛽^𝑖_11𝑝𝑝, 𝛽^𝑖_12𝑝𝑝, and 𝛽^𝑖_2𝑝𝑝 were directional between age groups. Results from applying the FLUed model using the Table 1 parameter settings for both influenzas are shown in Figures 8 and 9, respectively. Actual and simulated case data are shown in weekly units.

(a)

(b)

Figure 8 

Comparison of actual and simulated results for (a) weekly new infected cases and (b) cumulative cases normalized for 2008-2009 seasonal influenza A.

(a)

(b)

Figure 9 

Comparison of actual and simulated results for (a) weekly new infected cases and (b) cumulative cases normalized for 2009 swine-origin influenza A (H1N1).

Statistical correlation coefficient (CC) and coefficient efficiency (CE) test results of the two influenzas were 0.86 and 0.74 for the 2008-2009 seasonal influenza A and 0.77 and 0.36 for the 2009 swine-origin influenza A (H1N1), respectively. Figure 8 presents a plot of newly infected case fractions for the 2008-2009 seasonal influenza A in Taiwan between September 2008 and April 2009, normalized to total cases. Higher CC and CE values for 2008-2009 seasonal influenza A explain the similarities between the two epidemic curves. In Figure 9, we plotted fractions of new infected cases for the 2009 swine-origin influenza A (H1N1) from week 25 to week 52, also normalized to total cases. As shown in that figure, the number of actual cases decreased between weeks 37 and 48, followed by an increasing trend, thus resulting in a lower CE value. This dual-wave pattern is similar to global diffusion patterns associated with international travel [35–38]. The second wave in Figure 9 also coincides with Taiwan’s usual influenza season. In the temperate regions of the northern hemisphere, most influenza activity occurs from November through April. In their equivalents in the southern hemisphere, it occurs from April through October, while in the tropics all influenza viruses circulate at low levels, year round [39]. Thus, researchers pay special attention to travel patterns in the two hemispheres during their respective influenza seasons [40] since close contact among large numbers of individuals on airplanes and in airports puts travelers at higher risk of contracting influenza [39, 41]. Transmission is greatly accelerated when influenza spreads to the world’s 50 largest airports, which account for approximately 70% of all airline passengers [42]. Since we did not incorporate international travel, our FLUed model failed to capture the second wave; however, it did capture the peak period for the first (primary) wave.

Next, we evaluated the impacts of transmission rate reduction (Figure 10 and Table 2(a)) and the effects of various intervention policy scenarios on different reduced transmission rates (Figure 11) using the previously described model parameters. As shown in Figure 12, special emphases were placed on the peak numbers of infected cases (also called “new infected cases at epidemic curve peak”) and the time of peak infection (also called “week number of epidemic curve peak”). Public health officials have two goals: reducing the peak numbers of infected cases and delaying peak infection times. Both efforts are aimed at extending time for allocating and rearranging limited epidemic prevention and medical resources such as vaccines, antivirus drugs, medical personnel, and hospital negative pressure isolation beds, as well as for obtaining support from neighboring cities/countries to resolve issues tied to inadequate medical resources.

Table 2 

Observation index values according to different policy activation scenarios during swine-origin H1N1 influenza A outbreak in Taipei.

(a)

(b)

(c)

(d)

Figure 10 

Simulation results under various transmission rates. (a) New weekly infected cases at different transmission rates. (b) Cumulative new infected cases at different transmission rates. (c) Comparison of new infected cases at epidemic curve peak at different transmission rates. (d) Comparison of week number at epidemic curve peak at different transmission rates.

(a)

(b)

(c)

Figure 11 

(a) New infected cases at epidemic curve peak according to various intervention policy scenarios at different reduced transmission rates. (b) Total number of infected cases according to various intervention policy scenarios at different reduced transmission rates. (c) Week number of epidemic curve peaks according to various intervention policy scenarios at different reduced transmission rates.

Figure 12 

Dual-intervention policy evaluation index (peak numbers of infected cases and time of peak).

According to the weekly fraction data for new infected cases shown in Figure 10(a), epidemic scale and curve peak height were negatively affected by decreased transmission rate. As shown in Figure 10(b), the cumulative number of infection cases decreased nonlinearly as transmission rate decreased. Results of our comparisons of curve peaks for new infected cases at different transmission rates are shown in Figure 10(c); the peak for the original transmission rate shown in Table 1 is at the leftmost part of the graph. As shown, the transmission rate should be reduced by at least 50% to improve peak number suppression; an obvious decrease in peak number occurred at a 70% transmission rate reduction. Week numbers for epidemic curve peaks at different transmission rates are shown in Figure 10(d); note the nonlinear increases as transmission rate decreased—a positive and significant result for public health policymakers.

Data from various intervention scenario simulations are shown in Figures 11(a)–11(c). No differences in new infected cases at peak epidemic curve or total infected cases were observed for different intervention activation times (Figures 11(a) and 11(b)). However, epidemic curve peak week number was significantly delayed from 20 to 51 when intervention activation time was set prior to 50 cumulative cases with a minimum 70% reduction in transmission rate (Figure 11(c)). Activation time exerted a much weaker effect on peak timing at a 30% or lower transmission rate decrease. Our simulation results suggest that a reduced transmission rate between 50 and 90% was required for this intervention policy to have a positive effect on the 2009 swine-origin influenza A (H1N1) pandemic. It had little effect on a reduced transmission rate of 30% or lower, while incurring significant social costs. Our data indicate that intervention activation time did not result in a significantly reduced number of total infected cases or new infected cases at the epidemic curve peak but did have an obvious effect in terms of delayed peak time—a positive result for public health policy determination and preparation.

We compared differences among 2009 swine-origin influenza A (H1N1) initial outbreak areas and their effects on subsequent spreading (Figure 13, Tables 2 and 3). The Da’an district of Taipei was identified as a high-density hybrid area of businesses and residences, and the Wuri district of Taichung was identified as a low-density rural area. In the first (pre-swine-origin virus) scenario, the epidemic curve peaks appeared quickly when initial outbreaks occurred in Da’an rather than Wuri. The Taichung peak was significantly delayed when the transmission rate was reduced to 30%. Both locations had approximately the same number of new cases, but Da’an had a much larger number of total cases. After reducing the transmission rate from 50% to 30%, Wuri had a much later peak week, with no effect from intervention activation time on the total number of cases or newly infected cases in either location. These results suggest that less densely populated initial outbreak areas are more sensitive to intervention activation time—that is, the combination of early activation time and lower transmission rate leads to significant delays in epidemic curve peaks in less densely populated locations.

Table 3 

Observation index values according to different policy activation scenarios during swine-origin H1N1 influenza A outbreak in Taichung.

Figure 13 

Epidemic peak week numbers for Taipei urban and Taichung rural areas.

In addition to using peak infected case numbers and peak infection times, we also applied two quantitative indicators [31] to assess the performances of different intervention policies: prevention effect (6) and cost-efficacy (7). According to the first of these equations, the critical preventive effect indicator value is 1; we set this value at >1 to produce better effects since certain policies actually support the spread of an infectious disease at values of <1. We believe public health experts can identify optimal interventions by comparing different policies at this level. The cost-efficacy indicator in (7) was used to assess prevention effects per user-defined unit cost; its value was purposefully made positive to produce better effects (negative preventive effects occur at values below 0). Public health officials can use the same benchmark to make decisions regarding the best activation times for various intervention policies by concurrently comparing their preventive effects and cost-efficacy indicators. To activate an appropriate intervention policy at an optimal time, its preventive effect indicator must be as large as possible, and its cost-efficacy indicator must be as large as possible and positive. Consider=PreventioneffectPEpolicy𝑝Totalinfectedcaseswithoutactivatinganystrategy∈[],=×Totalinfectedcaseswith𝑝activated0,∞Cost-efficacyCEpolicy𝑝Totalinfectedcaseswithoutactivatinganypolicy−Totalinfectedcaseswith𝑝activatedTotalconsumedresourcecostsof𝑝⁻¹∈[−∞,+∞]. (6)(7)

After establishing prevention-effect and cost-efficacy definitions and usages, we used the two quantitative indicators and the parameters listed in Table 1 to compare five intervention policies at two initial outbreak areas (Da’an and Wuri) and on three policy application dates (weeks 1, 4, and 8 following the first confirmed diagnosis) for the 2008-2009 seasonal influenza A outbreak. The five intervention policies shown in Figure 14 were (a) vaccinating randomly chosen individuals (see the following Pseudocode 1); (b) identifying and inoculating individuals who had come into contact with infected individuals; (c) using public relations tools to strongly encourage hand washing and mask wearing by the general public; (d) enforcing home quarantines for infected individuals until they recover, including a minimum of eight days for individuals coming into contact with them; and (e) giving antiviral medicines to all individuals in advance of the expected epidemic.

(1) when policy is enacted do (2)   for each town 𝑖 in 𝑇𝑜𝑤𝑛𝑠 do (3)      for each age group 𝑝 in 𝐴𝑔𝑒_𝑔𝑟𝑜𝑢𝑝𝑠𝑖 do (4)       𝑆𝑖𝑝←𝑆𝑖𝑝×(1−𝐶𝑜𝑣𝑒𝑟𝑎𝑔𝑒_𝑟𝑎𝑡𝑒𝑖𝑝) (5)      end loop (6)   end loop (7) end

Pseudocode 1 

Vaccinating randomly chosen individuals (parameters list: Towns, Age_groups, Coverage_rate).

Figure 14 

Relationships of five intervention policies and the FLUed-related variables.

As shown in Figure 15 and Table 4, our simulation results indicate that the two most effective policies were (a) and (e), followed by (c). These policies were more effective when activated as quickly as possible; fewer and smaller differences in effect were noted when they were activated at later dates (1-2 months after the optimum date). Our main conclusions derived from the simulations were as follows: (a) the combination of hand washing and mask wearing by the general public was the most cost-effective policy and (b) using antiviral medicines in advance was more cost-effective than buying and using a mix of vaccines and anti-viral medicines.

Table 4 

Comparisons of (a) prevention effects and (b) cost-efficacies among five intervention policies.

(a)

(b)

Figure 15 

Comparisons of (a) prevention effects and (b) cost-efficacies among five intervention strategies.

4. Conclusion

Our goal in this paper was to integrate human commuting networks into multiple age-structured compartmental models to create a Four-layer Universal Epidemic Dynamics model (FLUed). The proposed model is capable of providing insights regarding the epidemic dynamics of emerging infectious diseases according to various interventions involving different initial outbreak locations, activation times, and policy suites. We believe that the FLUed model represents a convenient and effective tool for public health agencies responsible for initiating early responses to potential pandemics and for assessing transport intervention policies in outbreak locations. To build on this positive beginning, we plan to expand the four-layer framework to make it suitable for other acute diseases and to make it responsive to complex human contact structures.

Although our focus in this pilot study was on two influenza outbreaks in Taiwan, the FLUed framework is transferrable to other situations. The SLIR compartmental model in layer 1 represents a general epidemiological model for all droplet-transmitted respiratory infections. To simulate the dynamics of other infectious diseases, the SLIR compartmental model used in layer 1 can be replaced by other models such as SIS (Susceptible-Infectious-Susceptible) and SLIRS (Susceptible-Latent-Infectious-Removed-Susceptible). Furthermore, the age group and commuting interactions in layers 2 and 3 can be disassembled to meet the risk factor requirements of other emerging infectious diseases. As noted in an earlier section, the layer 4 network topology can be modified to meet the needs of different scales of link-node structures—for example, adding long-distance transportation networks to determine the impacts of shutting down railway or airline systems, or using different contact structures (e.g., mosquito-human) to model dengue fever, malaria, rabies, Japanese encephalitis, and other vector-borne and human-animal contact diseases. We believe that the FLUed model is suitable for detailed spreading scenario experiments and classroom teaching demonstrations of epidemiological issues associated with emerging infectious diseases.

Conflict of Interests

The authors declare no conflict of interests in the submission of this paper.

Authors’ Contribution

Chung-Yuan Huang and Tzai-Huang Wen contributed equally to this work.

Acknowledgment

This work was supported in part by a grant from the Republic of China National Science Council (NSC-101-2119-M-182-001 and NSC-101-2119-M-002-020).

Endnotes

1. A general-purpose epidemic simulation model must be able to (a) provide insights regarding the epidemic dynamics of emerging infectious diseases according to various interventions involving different initial outbreak locations, activation times, and policy suites; (b) assist with understanding the properties and efficacies of various intervention policies; (c) construct an effective, low-cost, and executable suite of intervention policies; and (d) reduce the difficulties and costs associated with learning epidemiological concepts [43].
2. An extendable epidemic simulation model must be able to alter the framework to make it suitable for other acute diseases and to make it responsive to complex human contact structures.
3. We used MathWorks MATLAB to create our proposed FLUed model. The use of Microsoft Excel to organize census and transportation data means that policy makers, public health professionals, and others who have less experience in specialized computer software will be able to generate various spreading scenarios with minimal assistance. For source code matching specific research requirements, please contact the corresponding author.
4. Nodes can represent any type of location or area chosen by a transportation researcher once the appropriate data are located and organized—for example, a region on an island, a town/city in a county, a county in a state, or a state in a nation.
5. In this paper, expert-based parameters are in a form allowing public health experts and computational epidemiologists to concentrate on estimating epidemic simulation parameters using influenza surveillance weekly report and related epidemic information provided by Taiwan Centers for Disease Control (TCDC).
6. As used here, epidemic velocity is the same as isopathic velocity, a term introduced by Berger and Luke [44] to quantify epidemic diseases. An isopath is a contour of constant disease level. As an epidemic progresses, isopaths expand outward from the initial infection site. Velocity is defined as the rate of isopathic movement in a particular direction, expressed as units of distance per time (e.g., miles/day).
7. Acceleration rate is defined as change in epidemic velocity over time.
8. In our proposed model, these assumptions can easily be changed using a GUI simulation console with specific input files.