The three basic elements of our national-level simulation model are (1) a previously developed stochastic agent-based model for disease spread at the community level; (2) detailed U.S. Census demographics and worker flow data for daily commuter traffic at short distances, and Bureau of Transportation Statistics data for the less frequent longrange travel behavior; and (3) high-performance parallel computing expertise in modeling millions to billions of particles on hundreds to thousands of processors. These three components, each of which we describe next in some detail, are brought together to provide a unique capability for a detailed modeling of disease spread in the United States population. (1) Community-Level Stochastic Simulation Model As the starting point for constructing our national simulation model, we employ a discrete-time, stochastic simulation model of disease spread within a structured 2,000-Germann, Kadau, Longini, and Macken Supporting Info: Page 2 of 37 2/26/06 person community. Similar models have been developed and applied previously to both influenza (1-4) and smallpox (5) . The model population is stochastically generated to match census-based nationwide distributions of age, household size, and employment status. Each person in the population belongs to one of five age groups: preschool-age children (0-4 years), school-age children (5-18 years), young adults (19-29 years), adults (30-64 years), and older adults (64+ years). Households consist of one to seven persons, with either one or two adults, and are grouped randomly into clusters of four households each, and further grouped into one of four non-overlapping neighborhoods, each containing approximately 500 people. Every person also belongs to a set of close and casual contact (also referred to as "mixing") groups, ranging from their household and household cluster (highest contact rates), to schools and workplaces, down to their neighborhood and the entire community (with the lowest contact rates, representing occasional interactions in malls, supermarkets, and churches, for instance). All preschool-age children are assigned to either a neighborhood daycare center, with 14 children on average, or to one of several smaller neighborhood playgroups, each with 4 children. Depending on their age, school-age children may belong to one of two elementary school groups (each shared between two neighborhoods, with 79 students each on average), to a community-wide middle school group (128 students on average), or to a community-wide high school group (average 155 students). These school contact groups are in general not actual schools, but rather representative of the typical daily interactions a student may have with classmates and other peers. According to US census data, 93% of children 5 -18 years old attend school, so we allow the remaining 7% to mix in the household, household cluster, neighborhood and community during the daytime. Working adults (restricted to those who are 19 -64 years old) belong to a work group of approximately 20 people. Although in reality many workplaces are larger than 20 people, we assume that workers make a contact of sufficient duration and/or closeness to transmit influenza virus with a subset of the entire workforce at that location. Transmission within each contact group is described by a contact probability c i (Table 3) , which may depend on the age of both the infectious and susceptible persons. This contact probability represents the likelihood (within each 12-hour period) of having a contact of sufficient duration and closeness for transmission of an infectious dose of Germann, Kadau, Longini, and Macken Supporting Info: Page 3 of 37 2/26/06 influenza virus to be possible between these two individuals in this social setting. The probability of transmission given such contact, P trans , is a single scalar number which multiplies each contact probability, allowing for a simple variation in contagiousness (typically represented by the basic reproductive number, R 0 ) without modifying the underlying social interaction network parameters. We do not allow for any seasonal or weekly variation in contact rates or transmission probability, and no births or non-flurelated deaths are included in our model. Each day, the probability of infection for each susceptible individual is computed based on the transmission probabilities for each potential infectious contact, p i = P trans × c i . If the infectious contact is receiving antiviral treatment, this transmission probability is further multiplied by (1 -AVE i ), where AVE i is the antiviral efficacy for infectiousness. Similarly, if they have been vaccinated, the vaccine efficacy for infectiousness VE i reduces the transmission probability by (1 -VE i ). The transmission probability p i can be further reduced for asymptomatic (yet infectious) contacts, as described in the next section. The probability of a susceptible person becoming infected is then computed as a product of all of the possible infectious contacts each day. Figure 3 illustrates this calculation for a susceptible adult (shown in blue) with one infectious child in the household (HH), one infectious workgroup (WG) contact, and three other infectious people in the wider community (Comm). The probability that this susceptible adult becomes infected is where c and a denote child and adult, respectively. (It is obviously trivial to generalize this to the case where the two community children have different levels of infectiousness due to therapeutic drugs or differing levels of severity, e.g. one symptomatic and the other subclinical.) A Bernoulli trial is conducted by generating a uniform [0,1] random number; if this number is less than P, the susceptible adult becomes infected and enters the latent phase of infection. If desired, the source of infection can then be determined by sampling from the relative contributions of each infectious contact to P (for instance, in this example an infection is most likely to be transmitted from the household child, but all 5 infectious contacts have a finite probability of being identified as the source). Germann, Kadau, Longini, and Macken Supporting Info: Page 4 of 37 2/26/06 We employ the same influenza natural history model as used previously (4), which for completeness is recapitulated in Fig. 3 . The main points are that the latent, incubation, and contagious period durations are each sampled from discrete distributions, with mean periods of 1.2, 1.9, and 4.1 days, respectively. (The contagious period includes both the slight difference between latent and incubation periods, as well as the standard postincubation period when symptoms appear in 67% of infected people.) Any infectiousness that is not accompanied by overt symptoms (namely, the post-latent part of the incubation period, if any, and the 33% of infected people who never develop symptoms) is assumed to be half as great as the infectiousness of symptomatic individuals, reducing the transmission probability p i by a factor of two. As before (1-5), we also allow people who become ill to withdraw from all contact groups except their household, with an age-dependent withdrawal probability and distribution of the number of days of illness before withdrawal for influenza taken from Elveback et al (1). The potential pandemic influenza strain was assumed to have an age-dependent attack rate pattern between the historical 1957-8 "Asian" influenza A (H2N2) (6) and 1968-9 "Hong Kong" influenza A (H3N2) (7) pandemic strains (see Table 4 ). For fitting purposes the attack rate pattern was calculated as an average of the final state of 500 independent communities that initially had 12 random infected people each (Fig. 4) . As a baseline the contact rates in households, small play groups, and large day care centers were taken from (3) where an H2N2 strain was modeled. However, since this attack rate pattern hits school-age children particularly hard (see Table 4 ), these rates had to be reduced by about a factor of 3 (this is also evident in (4), where a similar attack rate pattern was fit to a model specific for Thailand). The rest of the contact contribution was split between the remaining four contact groups (workgroups, household cluster, neighborhood, and community). Fine tuning to generate the contact probabilities shown in Table 3 was done by calculating the gradient vectors for the different age dependent attack rates with respect to the contact rate parameters, which gives a linear approximation of the dependence of the attack rates as a function of the contact rates. Germann, Kadau, Longini, and Macken Supporting Info: Page 5 of 37 2/26/06 Although the fitting was merely done for isolated communities, we find that the national model has a very similar attack rate pattern (see Table 4 ). The fundamental geographic unit in our model is the census tract, which is defined as a relatively stable geographic area with between 1,500 and 8,000 residents, with an optimum size of 4,000 people. In the 2000 , which is why a significant number of people (1.13 million, or 0.9% of the total workforce) were reported as working at locations 100 miles or more from their residence. We assume that such travel does not occur on a daily basis, and instead place these according to the age group-specific data in Table 5 . The average trip duration according to this data is 4.3 nights; we choose from a distribution between 0 and 11 nights according to the data in Table 5 . For the present implementation, each trip destination is a random neighborhood within a random community (including workgroup-only communities), which results in a simple "gravity" model with no distance information.  We assume that the pandemic influenza strain is introduced into the United States via arriving international passengers. Furthermore, we assume that by the time of this introduction there is an ongoing worldwide pandemic, so that there is no particular country or region which can be isolated (for instance, restricting arriving international flights from Southeast Asia). Consequently, we consider the 14 largest international airport gateways (14) in the continental United States (see Table 6 ), and introduce a small number of infected individuals each day. We do so by choosing a random tract and community within each county listed in Table 6 , and randomly infect between 0 and N individuals (chosen randomly from a uniform distribution) in that community. We take N proportional to the number of international arrivals at each airport, assuming 1-10 potential infecteds per 10,000 daily passengers. This represents a group of individuals, such as a family or business travelers, flying into the U.S. from the assumed pandemic Germann, Kadau, Longini, and Macken Supporting Info: Page 8 of 37 2/26/06 that is raging worldwide throughout the 180-day simulation. For most simulations we assume 2 infected persons per 10,000 daily international passengers, but the sensitivity to this choice and to other issues related to the seeding of infecteds is discussed below. The value of R 0 was calculated for different transmission probabilities by three different methods, yielding the results summarized in Table 7 . The first method was to average the number of secondary infections (omitting any tertiary infections) in 128,000 isolated communities that each had one random index case within the 2,000-person community population. (Using a smaller number of realizations, up to several thousand, led to statistical errors too great to determine R 0 within the desired ±0.1 precision.) An example of this calculated distribution and the resulting average R 0 are shown for P trans = 0.12 in Fig. 6 , and these results are denoted "random index case" in Table 7 . The second method was similar, except that separate R 0 first were calculated for index cases belonging to each of the 5 age groups (see Table 7 ). The overall R 0 was then calculated as an average of these age group-dependent R 0 values, weighted by the age dependent attack rate pattern for the respective transmission probability (referred to as the "attack rate pattern weighted index case" in Table 7 ). By doing this, the index case is more "typical" of those hit hardest by the outbreak, and as expected, this method slightly increases the value of R 0 (particularly for low R 0 ). The last method is an approximation based on the slope of the cumulative number of cases (15) and allows also for a time dependence of the reproductive number , where ν is the sum of the latent and infectious periods, which for our model is 1.2 + 4.1 = 5.3 days, f is the relative duration of the latent period (i.e., Fig. 7 ). Although there are large oscillations at early times due to the larger statistical errors (from fewer cases), it is clearly noticeable that R is largest at early times, and drops later. This is related to the fact that school children are particularly important spreaders in the initial stages of an influenza outbreak (see Fig. 7 ) due to their strong household and school interactions, which then naturally enhances R (see Table   Germann , Kadau, Longini, and Macken Supporting Info: Page 9 of 37 2/26/06 7). After about 30 days, the value for R stabilizes; in Table 7 we report (as "slope of cumulative number of cases") the average value between day 30 and the time when R begins to decline sharply because large parts of the nation are already affected. This approximation is in good agreement with the value obtained by the second method, involving only single-community, short-time simulations. We should note that the two former methods only can give an averaged static value as an estimator, whereas as the latter method can give information about the time development of R. Here, early-time fluctuations and enhancement of the reproductive number clearly demonstrate the difficulty in measuring this quantity from available data in a real epidemic. Furthermore, this behavior is even more complicated by the spatiotemporal spread of the epidemic, which causes local variations of R in time. It should not be too surprising that the averaged static value of R 0 for the community and national models are similar, since although the index case can interact with more than one community in the latter simulation, the number of effective interactions is the same, and in both cases the populations are completely susceptible. Upon activation of a TAP program, symptomatic individuals and their close contacts are treated with antiviral drugs, until a possibly limited national stockpile has been exhausted. We assume that X% of symptomatic cases can be identified, and that one day after the onset of illness, the sick individual is treated therapeutically and prophylaxis offered to his/her close contacts. Of these, we assume that 100% of household, household cluster, preshool, and playgroup contact are identified and treated, and that Y% of workgroup and elementary, middle, and high school contacts are identified and treated. For the present work, we will focus on two cases: X = Y = 60% or 80%, and refer to these as "60% TAP" and "80% TAP", respectively. As in other recent models of pandemic influenza (4), we use reported estimates of the antiviral efficacy for oseltamivir (16) (17) (18) (19)(20). Specifically, we assume that the antiviral efficacy for susceptibility AVE s = 0. course becomes ill, they complete their current course at the increased 2 tablets/day dosage. We assume that full antiviral efficacy is attained with the first tablet, and that there is no residual efficacy once the course has been completed. Two major uncertainties in modeling any vaccination program are how effective the vaccine will be (as even for endemic influenza, it is typically matched against a strain which is several months to a year old), and how quickly it can be produced, distributed, and result in an effective immune response. As yet, the efficacy of vaccination against a human-adapted avian strain is unknown. The immunogenicity of experimental vaccines has been measured; it has been found that a four-fold increase in antigen content above that of vaccines against human strains, and two vaccine doses, are required for a rise in antibody titre typically associated with protection (21). We assume that efficacy and immunogenicity are linearly related in our simulations. The second complication is related to timing. The time lag between vaccination and full effectiveness depends on the particular vaccine; for instance, a live attenuated vaccine may produce an antibody response within 1 day, while a killed vaccine may take two weeks. If multiple doses are required, the timescale can be significantly longer; for instance, two doses of a killed vaccine administered four weeks apart means that full efficacy may not be attained until six weeks after the initial dose. Rather than dealing with the specifics of any particular vaccine (including partial efficacy between the administration of the first and second doses), we simply combine this delay time with that for production and distribution, and refer only to the date at which vaccination becomes effective, which may be either before or after the outbreak begins. ** We consider two alternative distribution strategies, either randomly throughout the entire eligible population, or preferentially to children (with any remaining vaccine then distributed among adults). In either case, the eligible population consists of all individuals who have not been vaccinated and are not currently symptomatic. For ** This leads to a minor inconsistency, in that the eligible population is determined at the date at which full effectiveness is reached, and not at the earlier date of vaccination. However, this study is concerned with strategies to minimize the number of infected individuals, in which case the number of new symptomatic cases (who become ineligible for vaccination) between these two dates is negligible. Germann, Kadau, Longini, and Macken Supporting Info: Page 11 of 37 2/26/06 simplicity, we only consider two alternative production scenarios, either assuming the early distribution of a low-efficacy (e.g., a poorly matched) vaccine or the delayed production of a higher-efficacy vaccine. The well-matched vaccine is assumed to require two doses, and to have a vaccine efficacy for susceptibility VE s = 0.70 (with a reduced VE s = 0.50 for the elderly, age 65+), and a vaccine efficacy for infectiousness VE i = 0.80. The poorly matched vaccine has only VE s = 0.30 (for all age groups) and VE i = 0.50, and is assumed to require only a single dose (which would not be the case for an avian H5N1-based vaccine, for instance (21)). It is assumed that early production of the poorly matched vaccine allows for a vaccination program prior to the outbreak, resulting in a prior coverage of some fraction of the population (again, either uniformly or preferentially to children). For either vaccine, we assume a constant production and distribution rate of 4, 10, or 20 million doses per week nationwide, either starting as soon as two months before the first introduction, to as late as two months after the first introduction. The total production is also assumed to be limited to 50, 100, 250, or 400 million doses. Upon recognition of a pandemic strain in the U.S., one of the likely mitigation strategies is the closure of schools (22). We assume that this involves a total nationwide closure all of the school-related mixing groups in our model, and that this closure remains in effect for the duration of the pandemic. The affected mixing groups are the regular preschool-age playgroups, preschools, and elementary, middle, and high schools. All other contact rates remain unchanged. As a result of either a formal quarantine program, or voluntary changes in social and hygienic behavior in the event of a widespread pandemic, it is likely that normal contact behavior will be affected in times of crisis. Although this alteration is difficult to predict in advance, it almost surely will involve an increased tendency to remain at home rather than in large public places. To approximate this behavioral modification, we assume that the contact rates are cut in half for the community, neighborhood, work group, school, preschool, and playgroup mixing groups; household contact rates are doubled; and household cluster contact rates remain unchanged. As with the other mitigation Germann, Kadau, Longini, and Macken Supporting Info: Page 12 of 37 2/26/06 strategies, it is assumed that this alteration in normal behavior occurs nationwide, and lasts throughout the remainder of the epidemic. We consider reductions in both of our travel components: the daily workplace travel and irregular long-range travel. The first may be curtailed by voluntary increases in telecommuting, or in extreme cases by a nationwide work stoppage (with the exception of health care and emergency personnel, as described below). Reductions in longrange travel may also range from a component of the natural social distancing tendency, to an imposed quarantine or travel restriction program. We assume that long-range travel may be reduced to as little as 1% of the normal number of trips. In this section we explore the sensitivity to various components of the model, including how the epidemic is introduced, delays in implementing intervention strategies and the assumed effectiveness (or public compliance) with each intervention. Although each of these variations naturally leads to quantitative changes in the precise quantitative results, the basic conclusions presented in the main text are all relatively insensitive to any of the variations that we have explored. Since the model is inherently stochastic (including the mock population generation, the introduction of new infecteds, the daily disease transmission, and the intervention strategy components), in theory we need to run several realizations for each scenario, each with a different initial random seed. (For parallel runs, each processor uses a different seed to avoid having exactly identical communities anywhere in the simulation.) However, we have observed that between the large degree of spatial averaging over the 180,000 communities making up the national model, and the daily introduction of new infecteds which demands a robust intervention strategy, that there is almost no variation in the either the nationwide epidemic curve or final attack rate, but only subtle differences in the specific timing and geographic details of the epidemic spread. For instance, we compared eight different baseline (no intervention) realizations for R 0 = 1.6, differing only by the initial random number generator seeds. All eight simulations give Germann, Kadau, Longini, and Macken Supporting Info: Page 13 of 37 2/26/06 nearly identical final attack rates (four simulations give 32.62%, three give 32.63%, and one 32.64%), with the epidemic peaking between 115 and 120 days after the first U.S. introduction in all realizations. (2) Seeding of the epidemic In order to investigate the sensitivity against different introductions of infected individuals (size and spatial distribution of these individuals, as well as dynamic vs. static seeding), different amount of travel, as well as different random seeds (i.e. stochastic behavior), a series of runs with no intervention is presented for R 0 = 1.9 (similar results were obtained for other choices of R 0 ). Figure 8 shows the sensitivity to the dynamic seeding rate, from 1 to 4 infected persons per 10,000 international passengers each day at the major 14 US international air hubs, as well as a static seeding of 8 infecteds per 10,000 international passengers only at the beginning of the simulation (corresponding to 76 infected people introduced on day 0). It can be clearly seen that the size of the seed and the effect of static vs. dynamic seeding only shifts the epidemic curves, without affecting the shape or the overall attack rate (43.53% in all cases).
In many scenarios it is beneficial to trace the dynamics in a more fine-grained way especially during the initial or final stages of an outbreak, when person-to-person transmission processes dominate REF .
Mitigation strategies for pandemic influenza in the United
Proceedings of the National Academy of Sciences of the United States of America
Proceedings of the National Academy of Sciences of the United States of America
Biology Medicine