Introduction

In recent years, public health officials have been increasingly alert to the possibility of biological weapon use through an act of bioterrorism. Among potential pathogenic agents, smallpox poses one of the greatest risks. While smallpox has been eradicated in 19801, many experts believe that smallpox virus could exist outside the official institutes in the US and Russia. The smallpox virus, which is called variola major, is highly pathogenic and its intentional release into a population that is now largely susceptible (no one is certain about the residual protection of vaccination received over 30 years ago) can create a tremendous amount of social disruption and fatality. The possibility of a bioterrorist attack with smallpox is extremely limited and hampered by difficulties in using smallpox as a biological weapon; however, there is an obvious need for contingency planning and preparation. From a different perspective, smallpox is also used as the paradigmatic highly pathogenic agent that serves as a case study for the preparation of a bioterrorism event or a laboratory accident, as recently revamped by the engineered H5N1 debate2.

One of the main challenges in modeling intentional smallpox release is the choice of key epidemiological parameters and natural history of smallpox infection. The contagiousness of the disease in the present population, the release method of the virus and many options in public health response, such as ring and mass vaccination, mobility restrictions create a wide range of scenarios for which it is crucial to develop epidemic models able to gauge the actual threat of a smallpox bioterrorist attack. In recent years, several studies have focused on the study of smallpox transmission and control in urban settings and at the level of individual countries. Different modeling techniques3,4, ranging from compartmental models5,6,7,8,9,10,11,12 to highly detailed agent-based models13,14,15, generally agree in concluding that case isolation and vaccination, if timely implemented, will be sufficient to halt the ongoing transmission and contain the outbreak. All studies, however, consider the situation in which the public health system is effectively implementing the containment policies and is in possession of vaccine stockpiles. While these assumptions are likely to hold in many wealthy countries, it is hard to imagine the same coordination and availability of stockpiles worldwide. Furthermore, smallpox has long incubation and prodromal periods that can last more than two weeks. Finally, it is plausible that the correct diagnosis of smallpox cases will not be immediate, as doctors would initially not consider this eventuality. It is thus likely that two to four weeks could pass from the smallpox attack before a worldwide emergency is declared. This implies that during this time the disease might spread to other countries by means of traveling people exposed to the virus16.

The potential threat represented by the global effects of a targeted smallpox release event has yet to be analyzed with the use of explicit models. Here, building on previous work, we have developed a large-scale structured metapopulation model that accurately describes the worldwide spread of smallpox during the initial period of time between the occurrence of the intentional virus release and its detection by health officials. Thus we focus here on the level of worldwide diffusion of smallpox cases before any containment/mitigation policy can be implemented.

Smallpox's natural history is generally subdivided into three different stages: Latent, Prodromal and Rash1,17. Upon infection, individuals spend an extended period of at least Ld days incubating the infection6,10 in the Latent compartment, after which they progress with rate ε into the prodromal state12,13,15,16,18,19,20. The Prodromal period, with an average duration of γ−1 days6,9,12,15,18,19,20, is characterized by fever and flu-like symptoms before finally displaying the typical liquid filled blisters that characterize smallpox rash. Rash typically lasts Id days6,10 and is eventually followed by a late rash period15,20, after which recovery/death at a rate of μ occurs. The relative infectiousness associated with the prodromal period and the different rash stages changes considerably in the literature1,16,17. It is also plausible to consider that the transmission potential of the virus may change considerably because of the awareness of a smallpox outbreak. While infectiousness may be very high in the Rash stage, it is much simpler to detect and at least partially isolate those cases. Hospitalization of undetected cases at the beginning of the outbreak may lead to large transmissibility12,21,22, which is drastically reduced once the proper health infrastructure is set up. In the following, we consider only the earliest outbreak stages during which the public health system is unaware of the virus spreading. The basic reproductive number R0, gauging the transmissibility of the virus, is a parameter of the model and acts on the transmissibility of the various stages of the disease. We have considered, according to recent work12,13,15,23, that the Prodromal period contributes 10% of the disease transmissibility and that during the late rash period the infectiousness of individuals is reduced by 90% because the severity of the disease leads to isolation and identification of the cases. A key feature of our model is the mobility of individuals. It is therefore extremely important to associate the different stages of the disease with different mobility and travel capabilities. In particular, only a fraction pt of individuals in the Prodromal period are allowed to travel. Prodromal individuals capable of traveling are still allowed to travel with probability pr for Id days upon their progression to the Rash stage12. They are then withdrawn from the pool of mobile individuals13. The compartmental structure of the model, which includes the mobility classes, is reported in Fig. 1. In Table 1, we summarize the baseline parameters used in the study and the ranges considered for the sensitivity analysis of the model. All values and ranges considered are based on an extensive analysis of the literature and the parameters generally adopted in major simulation studies. The transmission dynamics is simulated with a multinomial stochastic transmission model and the individuals' mobility follows a stochastic dynamics too as detailed in the online supplementary information.

Table 1 Compartmental smallpox model parameters. Model parameters used in the baseline scenario and sensitivity analysis (in brackets) are reported
Figure 1
figure 1

Compartmental smallpox model.

Each susceptible individual in contact with an infectious case in Prodromal, Early Rash and Late Rash Stage, contracts the infection at rate rββ, β and rββ, respectively. Newly infected individuals enter a latency period during which they are not infectious yet and remain latent for a minimal duration of Ld days, after which they progress into the Prodromal Stage at per capita rate ε. Individuals in prodromal stage are divided between those who are able to travel, which occurs with probability pt and those who are restricted from traveling. Prodromal Stage is characterized by reduced infection transmissibility rββ and mean duration of γ−1 days. A proportion pr of individuals at the end of Prodromal Stage with traveling capability moves into Early Rash Stage to continue traveling while the rest – 1-pr – is withdrawn from traveling during this stage. Early Rash lasts for Id days and is characterized by the highest transmissibility β of the virus. After Early Rash, infectious individuals proceed to Late Rash Stage during which they are restricted from traveling and their transmissibility is reduced to rββ. This last stage of infection is followed by permanent recovery at rate μ.

Results

In order to model the spread of smallpox as occurring from an intentional release in the environment, it is crucial to imagine different release events19 and initial locations22. Similar to previous published studies, we consider two release scenarios. The first scenario depicts the intentional release from 5 individuals who successfully infect themselves with the virus in a major Western metropolis. These individuals move freely within the city (avoiding detection) and expose civilians to infection until they are removed from the population. The five individuals follow a deliberate release strategy and do not attempt international travel in order to limit the possibility of detection. As well, they do not seek medical care. In the model those agents follow just local mobility and keep on in transmitting the disease till they enter the late rush compartment. A second release event scenario deals with directly targeting the population. This can easily result in several successful infections if it were attempted in a context where tens of thousands of people are in close proximity. However, given the quantities of aerosol that would be required and the logistics necessary to distribute it to such large numbers while avoiding detection, we believe that this would be extremely unlikely and choose to focus on the more likely case of a small-scale initial infection as done in previous studies3,10,11,13. In particular we consider a case in which an aerosol version of the virus is dispersed in a closed environment, successfully exposing just 10 civilians. The victims are unaware that they have been exposed and continue to behave as they normally would. It is important to stress that the stochastic dynamics of the model naturally takes into account the generation of infections at different times and thus generates in each realization a varying number of exposed and infected people during the period preceding the detection of the outbreak by the health system. In both scenarios, we consider London, UK as the target city. In the supplementary information, we report results for other major locations in Europe and the US.

The conditions for the detection from the public health community of the intentional release are discussed in the methods sections. In particular, we assumed that a minimum number of 4 infected cases (sensitivity analysis between 2 and 8 cases) seek health care in the rush stage. Before detection the world is totally unaware of the intentional release, no containment/mitigation policy is therefore implemented. We consider the moment of detection as the moment at which, in the most optimistic scenario, the international community can start working on issuing and implementing the containment/mitigations contingency plans. This is a quite conservative assumption as medical doctors would hardly consider smallpox as a possible diagnosis until the onset of complications in the patients and therefore it is likely that the disease detection would require several more days to occur. Furthermore, the actual distribution of vaccine, implementation of travel restrictions, the organization of a coordinated effort of international contact tracing may require extra time. For this reason we provide here snapshot of the worldwide impact of the intentional release without considering mitigation/containment policies at the time of detection, as well as after 1 to 3 weeks after detection.

In Fig. 2 we report the number of countries affected and the smallpox cases observed outside the targeted country for the baseline transmission rate and for different delays after the initial outbreak detection. It is possible to observe that even at the moment of detection of the small-scale outbreak, with R0 = 5, the 50% reference range indicates that 2 to 4 countries have already imported at least one exposed individual. This is also observed with the number of exposed individuals dispersed in countries outside the target. The risk analysis for each country is composed of two components: the probability that a country will have at least one infected individual and the expected number of infected cases conditional to this event. In order to provide a global visualization of the risk worldwide, we provide risk maps at the level of single census area. In Fig. 3 we show a worldwide map that shows the probability of observing exposed individuals at the detection and after two weeks from detection of the release at the scale of the census area used in the numerical simulation. The statistics are obtained by analyzing 5,000 different stochastic realizations of the smallpox release with the same initial conditions. It is also important to stress that the reported outbreak probabilities refer to initial seeding events and small-scale outbreaks that may or may not lead to large-scale epidemics depending on local containment policies. In each different realization however the number of cases is generally localized in specific census areas. For instance if we consider Germany, we might observe a realization with 3 cases in Frankfurt, while in a different realization we may have 4 cases in Munich.

Figure 2
figure 2

Timeline of international smallpox spread.

Top) Intentional spread by five individuals avoiding detection and medical care. Bottom) Virus release with the exposure of 10 civilians. (A–B) Probability of observing an international outbreak (presence of cases outside the initial target country) at the time of initial detection and the following three weeks. (C–D) Number of affected countries conditional to the international outbreak. (E–F) Number of cases outside the target country conditional to global outbreak. Boxplots display the median, 50% and 95% reference ranges from 5,000 simulated events with the baseline parameters of Table 1.

Figure 3
figure 3

Worldwide location of risk areas in the case of the intentional release of smallpox virus in London.

The upper map reports the probability of observing exposed individuals at the time of detection and the lower map shows the same quantity two weeks from detection. The color code is in logarithmic scale and the results are at the finer scale allowed by the numerical simulations that depending on the region of the world ranges from 25 to 200 km. Maps were generated using ArcGIS.

According to our statistical analysis, not only is the risk never geographically restricted, even at the time of detection or just a week later when the first coordinated interventions might be feasible (see Figs. 2, 3), it also increases rapidly both in terms of exposure size and probability. The reason for such a rapid international spread is due to the fact that London serves as a major connecting point between the various continents, with daily connections to other parts of Europe, Africa, Asia and the Americas. However, other major hubs in the USA and Europe produce the same results, as shown in the supplementary information. At the early stage of the epidemic the worldwide exportation of cases from the target country can be assumed to be a Poisson process, thus assuming independence of each event and disregarding the possibility of cluster cases24,25,26. In this case the cumulative probability distribution P(T) to have a first positive event in any given country within time T, i.e., one infected individual traveling from the initially targeted country, will be

where λ(t) is the rate of positive case importation. The rate λ(t) depends on the size of the outbreak in the initially targeted country and the traveling rate of individuals, that is proportional to the passenger traffic between the two countries ω. In Figure 4 we show the local outbreak/importation of cases probability as a function of the daily passenger traffic of each country with the UK. At the time of detection the two quantities exhibit a high level of correlation (r>0.95, p<0.0001) and a roughly linear behavior. Indeed at the early stage of the outbreak the number of exposed cases in the initially targeted country is very small and the exponential term in Eq. (1) can be approximated as , yielding a probability of outbreak well approximated by a linear scaling with the passenger traffic: P(T)~ω. After three weeks from the detection it is possible to observe that the curve exhibits the exponential behavior of Eq. (1) as a function of traffic (figure not shown). This non-linear behavior is indeed at the origin of the inefficacy of travel restrictions in slowing considerably the global spread of epidemics. Even severe travel restrictions reducing traffic of 50% or more generally delay the importation of cases of only two to three weeks at the most as already pointed out in the literature24,27,28,29. However, in the case of a potentially catastrophic event such as a smallpox large-scale outbreak, drastic travel restrictions (to the extent of country isolation) may be appropriate in the framework of an international containment effort and deserve a separate careful study. This simple calculation shows clearly that the key parameter in the risk assessment of each country is the incoming traffic from the outbreak origin. A finer analysis considering catchment areas of specific airports and ground transportation and commuting is however required to achieve precise estimate especially within countries.

Figure 4
figure 4

Probability of observing smallpox cases as a function of the traffic from the targeted country.

The plot shows that initially the case importation/local outbreak probability in the top 20 countries at risk is well approximated by a linear scaling (solid line) with the traffic from the initially targeted country (high correlation coefficient r, with significant p-value). The size of the data points corresponds to the total population of the country. Data and scenarios are the same as reported in Table 2.

So far in this discussion, we have ignored any possibility of residual immunity from the global vaccination campaign that resulted in the eradication of smallpox 30 years ago. Several authors20,30 have pointed out that it is possible that as much as 20% of the global population remains immune, even though the long term effects of the vaccine are not well known. In order to test the effect of a residual immunity, we also inspect scenarios that allow for 20% of the population to be completely immune to contagion. Strikingly, we find that the overall results do not change significantly. There are roughly 40% less total cases, but they are still distributed across the globe, making a successful containment of the disease a major challenge. Further details and sensitivity analyses are presented in the supplementary information.

Many elements of uncertainty, some on the positive and some on the negative side, are present in the modeling assumption used in assessing the level of threat of the international spread resulting from the deliberate release of highly pathogenic virus. First of all, the transmissibility of the virus has an obvious impact on the international spread. More accurate modeling of refined population structure may enhance or hamper the international spreading depending on the mobility habits of the infected individuals. We do not include age structure or income differences in identifying travelers. The model does not include cluster events in the importation of infectious individuals such as those that may occur in the confined space of airplanes31. Similarly we do not include contact structure and index case setting that may lead to super-spreading event as observed for instance in the SARS epidemic32. It is clear, however, that time is a crucial factor. At the moment of the initial detection of an outbreak the number of infected individuals outside the target country is relatively small. A very effective contact tracing may lead to the timely identification and isolation of all cases before they can trigger large-scale outbreaks. On the other hand, one to two weeks of delay may lead to the impossibility of a 100% effective contact tracing of more than 50 individuals. The next step is, therefore, the study of the effect of quarantine and other mitigation and containment policies, especially in the event of outbreaks in less developed countries, at the international level.

Discussion

The present work indicates that a deliberate smallpox release is likely to assume an international dimension even before the epidemic is identified. We show through large-scale individual-based simulations that biological targeted attacks on a single city can result in the presence of exposed individuals in several countries before the health system is aware of the release and the ensuing outbreak. Some of the countries that could be affected may not have health infrastructures able to timely cope efficiently with the emergency dictated by a highly pathogenic virus outbreak. These findings highlight the presence of a systemic risk linked to modern human mobility that has to be faced by developing new contingency plans that consider the non-local dimension of epidemic spreading. The presented computational approach is a potential starting point for the definition of optimized containment schemes and mitigation policies that includes the international dimension of epidemics caused by highly pathogenic viruses. For reasons related to dual-use risk of the results contained in this paper, according to the comments of biosafety reviewers, we have removed quantitative data on risk probability and outbreak size in different scenarios. Those additional outputs can be shared with government officials and biosecurity researchers upon request.

Methods

GLobal epidemic and mobility model

We use the GLobal Epidemic and Mobility (GLEaM) model24,33,34,35,36,37,38, which is based on a high definition geographically structured metapopulation approach39,40. The model is composed of three layers. The first one, the population layer, integrates distinct census areas for a total of 3,362 subpopulations in 220 countries around the world. For each subpopulation, the number of individuals is obtained from “Gridded Population of the World,” a project by Columbia University that provides population estimates worldwide for cells of 15 × 15 minutes of arc41. The second layer of the model contains information about human mobility flows across the census areas. We consider both commuting flows, which is collected from various sources37 and airline traffic worldwide, which is provided by commercial databases42,43. Within each subpopulation, the infection dynamics takes into account the natural history and key parameters through a compartmental structure. The disease progression is simulated through a stochastic chain binomial and multinomial model44,45. Finally, the transportation layer allows for the detailed description of the mobility of exposed and infected individuals, thus allowing the stochastic simulation of the worldwide unfolding of the epidemic. Further details concerning the model are in the supplementary information and in Refs 35, 37, 38. As we work with a fully stochastic model, we have considered, for each considered scenario and parameters' set, a total of 5,000 stochastic realizations. The subpopulation structure covers explicitly almost 6 billion individuals with a time resolution of one day. Each set of realizations produces 27 GigaBytes of data and runs on 256 CPU cores.

Detection

What we are interested in is the quantitative assessment of the risk of the internationalization of the epidemic and the number of countries involved. Clearly, this risk assessment depends on the time elapsed from the biological attack. The point of interest is therefore the time at which the international community is able to issue a worldwide alert and start implementing containment and mitigation policies. This includes the ability of an effective contact tracing, the deployment of vaccine stockpiles for ring vaccination, travel restrictions, etc. Theoretically, the earliest time at which the detection can occur is when the first person in the rash stage is correctly diagnosed. Due to the current rarity of the disease and the difficult diagnostic, it is likely that an alert will only be possible after several civilian cases have already occurred. We assume that 4 individuals are required to enter the rash stage before successful detection is feasible and an alert issued46,47,48,49 and make a sensitivity analysis in the range of 2 to 8 individuals. We are, therefore, measuring the distribution of the number of countries affected and the specific risk posed for each individual country at detection time. In addition, from the successful detection of the smallpox outbreak to the implementation of effective contact tracing and worldwide coordinated response, including the deployment of a vaccine, most of the experts consider a time window ranging from 1 to 4 weeks. For this reason, we report data also for the three consecutive weeks (four consecutive weeks in the Supplementary Information) following the initial outbreak detection.