Modeling the Worldwide Spread of Pandemic Influenza: Baseline Case and Containment Interventions

Background The highly pathogenic H5N1 avian influenza virus, which is now widespread in Southeast Asia and which diffused recently in some areas of the Balkans region and Western Europe, has raised a public alert toward the potential occurrence of a new severe influenza pandemic. Here we study the worldwide spread of a pandemic and its possible containment at a global level taking into account all available information on air travel. Methods and Findings We studied a metapopulation stochastic epidemic model on a global scale that considers airline travel flow data among urban areas. We provided a temporal and spatial evolution of the pandemic with a sensitivity analysis of different levels of infectiousness of the virus and initial outbreak conditions (both geographical and seasonal). For each spreading scenario we provided the timeline and the geographical impact of the pandemic in 3,100 urban areas, located in 220 different countries. We compared the baseline cases with different containment strategies, including travel restrictions and the therapeutic use of antiviral (AV) drugs. We investigated the effect of the use of AV drugs in the event that therapeutic protocols can be carried out with maximal coverage for the populations in all countries. In view of the wide diversity of AV stockpiles in different regions of the world, we also studied scenarios in which only a limited number of countries are prepared (i.e., have considerable AV supplies). In particular, we compared different plans in which, on the one hand, only prepared and wealthy countries benefit from large AV resources, with, on the other hand, cooperative containment scenarios in which countries with large AV stockpiles make a small portion of their supplies available worldwide. Conclusions We show that the inclusion of air transportation is crucial in the assessment of the occurrence probability of global outbreaks. The large-scale therapeutic usage of AV drugs in all hit countries would be able to mitigate a pandemic effect with a reproductive rate as high as 1.9 during the first year; with AV supply use sufficient to treat approximately 2% to 6% of the population, in conjunction with efficient case detection and timely drug distribution. For highly contagious viruses (i.e., a reproductive rate as high as 2.3), even the unrealistic use of supplies corresponding to the treatment of approximately 20% of the population leaves 30%–50% of the population infected. In the case of limited AV supplies and pandemics with a reproductive rate as high as 1.9, we demonstrate that the more cooperative the strategy, the more effective are the containment results in all regions of the world, including those countries that made part of their resources available for global use.

city we assume a standard compartmentalization of the population in which every individual belongs to a certain compartment; the compartments used allow for a detailed description of a typical evolution of the actual disease, by including the incubation period and the asymptomatic course of influenza [7,8]. Individuals in each city can change compartment by following the infection dynamics and are also allowed to travel from one city to another by following the routes of the airline transportation network.
Details of the compartmentalization, in the baseline case and in the intervention scenario, are shown in the diagrams of Figure 2 of the manuscript.
In each city j the population is given by N j = S j (t) + L j + I t j (t) + I nt j (t) + I a j (t) + R j (t), where S j (t), L j (t), I t j (t), I nt j (t), I a j (t) and R j (t) represent at time t the number of susceptible, latent, symptomatic traveling, symptomatic not traveling, asymptomatic infected and recovered individuals, respectively. The model in the baseline scenario is described by the following set of discretized epidemic Langevin equations: S j (t + ∆t) − S j (t) = −β (I t j + I nt j + r β I a j )S j N j ∆t + β I t j S j N j ∆t η βt,j (t) + + β I nt j S j N j ∆t η βnt,j (t) + r β β I a j S j N j ∆t η r β ,j (t) + Ω j ({S}) (S1) L j (t + ∆t) − L j (t) = β (I t j + I nt j + r β I a j )S j N j ∆t − L j ∆t − β I t j S j N j ∆t η βt,j (t) + − β I nt j S j N j ∆t η βnt,j (t) − r β β I a j S j N j ∆t η r β ,j (t) + + µI a j ∆t η µa,j (t) + Ω j ({I a }) (S5) R j (t + ∆t) − R j (t) = µ(I t j + I nt j + I a j )∆t − µI t j ∆t η µt,j (t) − µI nt j ∆t η µnt,j (t) + − µI a j ∆t η µa,j (t) + Ω j ({R}), where η βt,j , η βnt,j , η r β ,j , η ,j , η µt,j , η µnt,j , η µa,j are statistically independent Gaussian random variables with zero mean and unit variance and ∆t represents the time step set to 1 day.
The set of Langevin equations are coupled by the stochastic transport operator Ω j which describes the movements of individuals traveling from one city to another. Following ref. [1], the transport term is a function of the air traffic flows between two connected cities and represents the net balance of passengers who left and entered city j.
The number of passengers of category X traveling from a city j to a city is an integer random variable, in that each of the X j potential travellers has a probability p j = w j /N j to go from j to . In each city j the number of passengers ξ j traveling on each connection j → at time t define a set of stochastic variables which follows the multinomial distribution where (1 − p j ) is the probability of not traveling, and (X j − ξ j ) identifies the number of non traveling individuals of category X. We use standard numerical subroutines to generate random numbers of travellers following these distributions.
The transport operator in each city j is therefore written as where the mean and variance of the stochastic variables are ξ j (X j ) = p j X j and Var(ξ j (X j )) = p j (1 − p j )X j . Direct flights as well as connecting flights up to two-legs flights can be considered.
In addition to the previous fluctuations, the transport operator is in general affected by fluctuations coming from the fact that the occupancy rate of the airplanes is not 100%. To take into account such fluctuations, we assume that on each connection (j, ) the flux of passengers at time t is given by a stochastic variablẽ where α = 0.7 corresponds to the average occupancy rate of order of 70% provided by IATA and η is a random number drawn uniformly in the interval [−1, 1] at each time step.
This system of differential equations can be numerically integrated to obtain the evolution of the infection dynamics inside each city. A major issue in the integration of the epidemic equations is the consideration of the discrete nature of the individuals and the preservation of the symmetry of the noise terms. For this reason we adopt an integration procedure that separates each variable X j into its integer part [X j ] and its non-integer remaining partX j and we write the corresponding time evolution for these variables [19,1]. This procedure is implemented at all time steps in each urban area j and provides a discrete output which preserves the correct integration of the equations. In addition, it is worth noticing that the derivation of the Langevin equations from the master equations describing the evolution of the probabilities associated to the infection dynamics is valid under the assumption of large population in each compartment for each city [18,16,17]. Therefore for the sake of mathematical rigor, in each compartment with less than 10 3 individuals we describe the evolution of the disease by fully considering discrete binomial or multinomial processes. This condition is checked at all time steps for all compartments in all urban areas. In particular two different kind of processes in the infection dynamics have to be considered: the generation of new cases (i.e. the generation of latents through the contact of susceptible and infectious individuals) and the transition of individuals from one compartment to another (e.g. from latent to infectious, from infectious to recovered).
For the first class of process it is assumed that each susceptible will become in contact with an infectious individual with probability βI j (t)/N j , so that the number of new infections in city j is extracted from a binomial distribution with probability βI j (t)/N j and number of trials S j (t).
Analogously, the number of individuals changing compartment in a transition process (e.g. in city j: L → I with rate ) is extracted from a binomial distribution with probability given by the rate of transition (in the previous example, ) and number of trials given by the number of individuals in the compartment at time t (in the previous example, L j (t)).
In case of multiple generation of new cases (due e.g. to susceptible individuals coming into contact with asymptomatic and symptomatic infectious individuals) and transition processes from a single compartment to more than one (e.g. latent becoming either symptomatic or asymptomatic), the binomial is substituted by a multinomial distribution, whose parameters are given by the probabilities associated to the single processes considered and the population of the given compartment. As an example, the number of latent becoming either asymptomatic or symptomatic traveling/non-traveling is extracted from a multinomial distribution based on the out of the total pool of latent individuals at time t, L(t).
When the assumption of large populations is satisfied, the epidemic evolution is then obtained through the integration of the Langevin equations, as described above. Results are in agreement with those obtained from a fully discretized version of the model in which all processes of the infection dynamics are described through multinomial distributions.

Infectious period under AV treatment
The efficacy of the AV treatment is modeled through a reduction of the infectiousness of treated individuals and a reduction of the total infectious period [7,8,10,9,20,22]. Since the compartmentalization considered does not allow to track the time since infection occurred for each individual, we consider a reduction of x days of the average infectious period. We assume that infectious individuals are allowed to recover before eventually being detected and administered an AV dose. Therefore, under AV treatment the transition from infectious individual to infectious under treatment -I t,nt to I AV -is assumed to occur with rate (1 − µ)p AV , while the rate of transition directly to the recovered compartment is left unchanged. The average period µ −1 AV spent in the I AV compartment is then obtained by reducing the average total infectious time by x∆t: where the first term in the l.h.s. of the equation represents the average time spent by a treated individual in the I t,nt compartment. From the equation above we obtain the following expression for the average time spent in the treated compartment: which clearly depends on the distribution rate p AV . Indeed, a large value of p AV refers to a prompt identification of cases and subsequent drug administration. Plugging into eq. S11 the values of the parameters considered (µ), we obtain different values of the rate of transition µ AV for different reductions x and rate distributions p AV .
Some values of x and p AV lead to a value of µ −1 AV smaller than 1. Since our time scale is ∆t = 1 day, in those cases we fix µ −1 AV to its smallest possible value, i.e. µ −1 AV = 1, and reduce the infectiousness β(1 − AV E I ) of the treated individuals by an appropriate scaling factor s β in order to effectively compensate the longer average period spent under treatment, imposed by time scale constraints. In this way, it is also possible to study reductions x smaller than 1 day (see in the following the sensitivity analysis on the reduction of the infectious period under treatment).

Estimate of the reproductive ratio with AV treatment
While it is possible to give an analytic expression of the reproductive rate for the compartmentalization adopted in the baseline case as provided in the manuscript, the same task becomes more difficult when AV intervention is considered because of the time delays due to the interplay of travel and AV use. In the following we provide a rough estimate of the reproductive rate under the assumption that AV intervention is put in place since the start of the pandemic (quite unrealistic), with no additional delays. In this case, the effective R 0,AV would be: where the r.h.s. of the equation represents the contributions of (i) asymptomatic individuals, (ii) symptomatic non-treated individuals and (iii) symptomatic treated individuals, respectively. In order to have a quantitative estimate of the value of R 0,AV with respect to the reproductive rate R 0 in absence of AV, we can compute the ratio: (S13) For a distribution rate of p AV = 0.7/day and assuming a reduction of the infectious period x = 1 day, the value of the ratio is R 0,AV /R 0 0.77, so that for R 0 = 1.5 the estimate for the effective R 0,AV gives a value 1.15, i.e. slightly larger than 1. However, this estimation does not take into account the fact that at initial stages of the original outbreak and of all secondary outbreaks the reproductive rate is not lowered by the use of AV treatment since the intervention takes place with a characteristic delay. The interplay of the travel and of the delays with which intervention strategies are put in place allow therefore the virus to spread out of the source, even in case of a high distribution rate of AV doses -p AV = 0.7 -and reproductive rates R 0 1.5.

Initial conditions and geographical subdivision
At the beginning of the simulation, the entire population is considered to be susceptible. The population vulnerability to an H5N1-like pandemic is indeed expected to be universal, since the human immune system will have no pre-existing immunity to a new virus [11]. We simulate the start of the pandemic influenza on several different initial dates, in order to assess the impact of

Seasonality
In order to take into account the high influenza incidence during winter period and low incidence during summer season, we include seasonality in the infection transmission of the model. We follow the approach described in refs. [13,14,15] by using a transmission parameter β which varies in time and depends on the geographic zone. The world is divided into three different climate zones -Northern hemisphere, Tropical zone and Southern hemisphere -as outlined by the Tropic of Cancer and the Tropic of Capricorn (see Figure S1).
To each city, we assign a scaling factor for the transmission parameter which depends both Figure S1: Subdivision of the world into three climate zones -Northern hemisphere, Tropical zone, Southern hemisphere -as outlined by the Tropic of Cancer and the Tropic of Capricorn.
on the time of the year and on the climate zone to which the city belongs to (see Table S2); its value models the seasonal variations of influenza. Following [15], we assume an enhancement of 10% of the transmission parameter β during the influenza season (November-February in the Northern hemisphere and May-August in the Southern hemisphere) and a reduction to as low as 10% during the non-influenza season (June-July in the Northern hemisphere and December-January in the Southern hemisphere). A linear interpolation between the two values is assumed to compute the scaling factor value for the remaining months. For all cities in the tropical climate zone, we assign a scaling factor equal to 1.

Stochastic noise
In Figure S4 we report the same results shown in Figure

Sensitivity analysis
In the following we show the results obtained from the sensitivity analysis with respect to different factors characterizing the global spread of the pandemic as illustrated in the manuscript.
More precisely, we study the effect of different temporal and geographical settings for the start of the pandemic in both the baseline and the intervention scenario; we also examine the role that the AV protocol assumed, the value of the reproductive ratio and the time delay of the intervention have on the outcoming results.  Analogously to the results obtained for the effect of the AV protocol, also in this case we observe an increase in the global attack rates after one year, with probability distributions for the number of infected countries almost unchanged. However, the increase observed is smaller than the one obtained by changing the distribution rate p AV . Figures S7 and S8 show that a large delay in the intervention could be overcome with a prompt detection of symptomatic individuals and a high rate distribution of doses.

Intervention with AV therapeutic treatment: limited AV supplies
This section will analyze the effect of several different factors characterizing the intervention strategy through the use of limited AV supplies in the general scenario. The AV protocol assumed is p AV = 0.7/day.

Starting date
If the pandemic originates in Hanoi in April, no striking differences are observed with respect to the case discussed in the manuscript, where the starting date is October. As Figure S9 shows, we still observe a strong reduction of the number of cases per 1000 persons and a delay in the prevalence profile when the intervention is in place with respect to the baseline curve.

Initial seed
Results obtained for a pandemic seeded in Bucharest in October ( Figure S10) are in agreement with the findings presented in the manuscript. In particular, it is possible to show that

Redistribution of shared AV resources
Here we analyze a different method of redistribution of the shared AV resources provided by the prepared countries. Instead of being used upon need, AV stockpiles are preemptively distributed to each of the unprepared country proportionally to its population. Figure S11 shows that this solution is less effective with respect to the distribution of shared resources from a unique global stockpile for the country use when in need. Indeed, while some

Reduction of the infectious period during AV treatment
Results presented in the manuscript assume that the reduction of the average infectious period induced by the AV treatment is of 1 day. Here we study the impact of the efficacy of AV drugs as measured by its reduction of the infectious period. Following the details of the infectious period distribution as outlined in a previous section, we assume a reduction of x = 0.5 days.
While the global attack rates are larger than the ones obtained with x = 1 day, as expected, cooperative strategies are still more efficient than the uncooperative one in reducing the number cooperative strategy II), are not able to consistently reduce the attack rate as in the case when p AV = 0.7/day, due to the high traffic air connections with South East Asia which is experiencing a large outbreak. This clearly shows that, given the same amount of AV doses along with the same distribution of resources, it is preferable to treat the largest possible percentage of ill individuals in the unit time. Since we are dealing with finite and limited resources, one could actually expect that a better protocol would be to distribute a lower percentage of doses, in order to keep the available supplies for as long as possible. Results show that this is not the case. Indeed, if the intervention assumes p AV = 0.5/day instead of a higher rate value, an explosion of cases will occur in some regions so that the resources will not last longer and the observed attack rate will be higher. Finally, if the reproductive rate is higher, the most important advantage of the ccoperative strategies is in the delay of the epidemic peak. Indeed, in this case many regions would experience the peak within the first 12 months, so that the global attack rates are not strikingly reduced as observed for a lower R 0 , but the vaccine development would still benefit of an additional delay of 1-2 months with respect to the uncooperative strategy.       Figure S16: Effect of different time delays of intervention in the countries experiencing secondary outbreaks. Here we consider a pandemic with R 0 = 1.7 starting in Hanoi in October with an assumed AV protocol of p AV = 0.5/day. From top to bottom, results for a delay of 7, 14 and 28 days are shown.

Travel limitations
While the enforcement of travel restrictions to and from infected areas is considered not applicable in most situations [11], it is reasonable to imagine that the air traffic volume in case of a pandemic will be reduced either due to the voluntary choice of individuals or for the applications of non-medical intervention measures. In Figure S17 we show what would happen in case the number of passengers on each connection is reduced by 20% or 50%.   Figure S17: Travel limitations as a non-medical intervention measure for a pandemic with R 0 = 1.9 starting in Hanoi in October Although the traffic limitations imply a considerable reduction of air travel, the two intervention measures do not appear to be efficient neither in reducing the epidemic impact, nor in delaying the prevalence peak, in agreement with the results of refs. [20,21]. The reduction in the air traffic is not able to contain the pandemic at the source. Indeed, as the outbreak occurs in the seeded country a larger number of cases is generated, so that the expected number of people traveling out of the country becomes considerably large though the probabilty of traveling has been drastically reduced. As soon as some infectious individuals travel from the source of the infection to another country, they will seed other regions with the infection, thus starting multiple outbreaks in different geographic zones (see the map reported in Figure 6 of the manuscript as an example). At that stage, the travel limitations are no longer able to considerably slow down the pandemic evolution.

Supporting Videos
Videos S1 and S2 represents the average simulated baseline scenarios corresponding to R 0 = 1.7 and R 0 = 2.3, respectively. The pandemic starts in Hanoi (Vietnam) in October. The color code represents the average number of cases per 1000 in each country on a logarithmic scale, from 10 −2 (green) to 1000 (red). The videos show a time period of two years (video S1) and 14 months (video S2) since the start of the pandemic. Every snapshot represents the situation at the beginning of each month.
Video S3 reports the results obtained for the first 12 months spread of a pandemic with R 0 = 1.7 in case the uncooperative strategy is applied, where only few countries have AV stockpiles available for 10% of their population. Initial conditions are the same as in videos S1 and S2.
Video S4 represents the cooperative scenario II described in the main text, where 1/5 of prepared countries' stockpiles is globally redistributed. Initial conditions and infectiousness of the virus are the same as in the previous videos. The video shows the evolution during the first 12 months from the start of the pandemic.