Fractional diffusion emulates a human mobility network during a simulated disease outbreak

From footpaths to flight routes, human mobility networks facilitate the spread of communicable diseases. Control and elimination efforts depend on characterizing these networks in terms of connections and flux rates of individuals between contact nodes. In some cases, transport can be parameterized with gravity-type models or approximated by a diffusive random walk. As a alternative, we have isolated intranational commercial air traffic as a case study for the utility of non-diffusive, heavy-tailed transport models. We implemented new stochastic simulations of a prototypical influenza-like infection, focusing on the dense, highly-connected United States air travel network. We show that mobility on this network can be described mainly by a power law, in agreement with previous studies. Remarkably, we find that the global evolution of an outbreak on this network is accurately reproduced by a two-parameter space-fractional diffusion equation, such that those parameters are determined by the air travel network.


Introduction
Characterization of the spatial spread of a disease is an essential problem for disease control and elimination. Epidemiological models that incorporate spatial dispersion have been used for numerous diseases such as measles [1,2], dengue [3,4,5], seasonal influenza [6], and pandemic flu [7,8,9]. A recent review of the topic [10] identified five opportunities for including spatial effects in epidemic models, including the issue of long-distance transport. The underlying spatial process for many epidemiological models is often assumed to be classical diffusion, gravity, or radiation flux. However, these are not necessarily valid for human mobility networks with long-range movements, as was shown for the circulation of US bank notes by Brockmann et al. [11].
Their seminal paper used ambivalent fractional-diffusion as a generalized model of transport processes with extreme events, which matched observed trajectories of bank notes. A similar approach, space-fractional diffusion, arises from a generalization of a random walk when the distribution of step sizes in the walk is described by a power law, x −(1+α) with α ≤ 2, rather than a Gaussian function of the step sizes. This type of model has general application across fields from plasma physics [12] to ecology [13,14]. Mathematically, this heavy-tailed distribution of step sizes with an infinite variance produces superdiffusive dispersion, also known as a Lévy flight. Waiting times between steps are taken to be λ-scaled within a Poisson process. In practice, a natural process can only be roughly approximated by a Lévy flight. In some applications, a truncation of the Lévy flight or a similar Lévy walk can be a good approximation [15,16,12]. Such models cannot expected to account for all the detailed connections in a non-local network, but may provide a useful and efficient macroscale approximation for a practically heavy-tailed (long-range) distribution of step sizes.
Human displacements occur on many scales through channels ranging from walking to driving to flying. The relative importance of these channels for spatial epidemic models depends on the density and territorial extent of a society, as well as cultural context and socioeconomic development. We decided to focus on the scheduled air traffic network that determines longrange transport in developed countries with large land area, such as the United States, China and India. Clearly, long-range displacements in many contexts will be important for understanding how diseases spread spatially [17]. We took an empirical approach to measuring air travel networks, rather than attempting an analytical result, as other have before [18]. These air travel displacements are instrumental in the fast propagation of a disease outbreak across long distances in certain societies. They are a major source of non-diffusive behavior as long flights bestow heavy-tailed distributions of step-sizes, transporting infected individuals long distances effectively instantaneously.
For the current work, we ran and analyzed spatial simulations of disease dynamics using a flight schedule database to determine passenger transition rates between areas serviced by commercial flights. Each airport acts as a catchment area for a population subject to a susceptible-infected-susceptible (SIS) disease dynamics. We limited our mobility model to air travel since we are interested in highlighting the long-range behavior that is relatively insensitive to local commuting [19,20]. Others [21] have already performed stochastic simulations of a realistic pandemic using global mobility data. Our main intention here was to use the dense air traffic network in the continental United States as an example of a human mobility network that can be approximated with a simple space-fractional diffusion model for transport. Others have recently examined a variety of network structures for disease dynamics [22,23]. Additionally, we extended the stochastic fractional-diffusion work presented in [24] with a more accurate, spatially second-order fractionaldiffusion method. Our method is related to fractional reaction-diffusion systems previously explored analytically and numerically [25,26,27,28]. Our results can feasibly be applied by epidemiologists for other transport networks that have heavy tails in the size of displacements between nodes [29].
Our presentation is organized as follows. First we show our analysis of commercial air travel on several countries' national flight networks before focusing on the USA. We extracted a superdiffusive transport exponent from the distribution of step sizes in the network data. We then describe our version of a stochastic SIS simulation to model seasonal influenza, using USA air travel data to form a weighted adjacency matrix connecting populated areas exposed to the disease. We found these results to be in agreement with known seasonal SIS behaviors, validating our new simulations. Next we propose a one-dimensional, space-fractional-diffusion equation as a macroscale substitute for simulating the United States' commercial air traffic mixing network. Remarkably, we found that our space-fractional approximation reproduces the national-average SIS infection curve best when the fractional exponent is similar to the superdiffusive exponent measured from the national flight network. In cases such as the one studied here, where long-range displacements form an important component of mobility, fractional diffusion is a promising macroscale approximation for general measures of epidemics.

Determining random walk parameters from a flight database
In order to approximate a complex human mobility network with a simple random walk model, the statistics of displacements on the network must be determined. We used a commercial air travel database (OAG Aviation Worldwide Limited) to find the probability density of displacements, in miles, for individual travelers. We set the scale of waiting times between steps at one week by aggregating air travel weekly. This should be a reasonable approximation to the travel habits of the average traveler. When aggregating the full set of flight data for several countries' air travel networks, a power law distribution of displacements is readily observed. Countries with large land areas and a dense airport network tend to have a heavy-tailed distribution of air travel step sizes, as shown in Figure 1(a). Here we have accounted for the number of seats in each aircraft and assumed all flights are full. This implies that we are taking an upper bound for mobility by air travel. Other authors have conducted more detailed analyses of air travel networks [30], but our focus was on reproducing large-scale epidemic dynamics with a twoparameter random walk model based on the distribution of step sizes.
Relative to those from the USA and China, displacement distributions from Brazil and India have a shallow tail with a shorter flight range due to the geographic area and distribution of population centers in these countries. The USA network step-size distribution has a superdiffusive power-law exponent 1 + α ∼ 1.5 ± 0.2 for a significant fraction of displacements in the main part of the tail of the distribution. We fit this power law within the range shown by the black lines using the Nelder-Mead method through the open Matlab toolbox Ezyfit. This power law is truncated to an exponential decay for steps larger than 1500 miles for the USA. We note that a superdiffusive scaling with α = 0.6 was also inferred from the dispersion of US bank notes [11]. India, China and the USA show similarities in the shape of the step-size distribution up to the scale of the cutoff. Since we were interested in reproducing this human mobility network with a very simple model, we approximated the US displacement distribution with a superdiffusive power law. For simplicity, we chose to neglect the effects of the naturally occurring truncation that have been discussed for other applications [16,31].
There are two parameters necessary for the partial power law fit indicated in Figure 1 Here Γ is the extension of the factorial function, π is the transcendental number, t = 1 since our random walk statistics are time-invariant. This formula [32] allowed us to estimate the order of magnitude of D α , which we found to be D α ∼ 1 when α ∼ 0.4 and D α ∼ 10 when α ∼ 0.7. This well-grounded, if inexact, calculation of D α is supported by our parameter scan in the final results section. Further discussion of the fractional diffusion equation is given in the Methods. The distributions of airport connectivity for the same four countries are given in Figure 1(b), where connectivity for a given node is simply the number of other nodes with a connecting flight to that node. This plot makes clear that the United States' network is roughly ten times denser than the next example, China, though the shapes of the curves are similar. With 482 nodes, the US air travel network is significantly denser than any other single nation examined. This network was therefore the focus of our study to find a plausible fractional diffusion approximation for a human mobility network superimposed with disease kinetics.

Basic model for influenza outbreaks on an air travel network
As we examined the macroscale dynamics of disease transport, we required a prototypical epidemiological model. We used the susceptible-infectedsusceptible (SIS) differential equations (see Methods) with β being the rate of infection in units of population per week, and γ the rate of recovery in the same units. In this standard model, the basic reproduction number is R 0 = β γ S 0 N , where S 0 is the initial number of infected and N = S + I is the total population size, including infected individuals numbering I (see Methods). Since we considered an outbreak at a single location, the fraction S 0 /N for the entire network is effectively unity. While our model was sufficiently generalized for any disease with SIS behavior, we made considerations for influenza and influenza-like illnesses (ILI). Since the flu virus evolves from season to season to reinfect populations, we treated the dynamics as effectively SIS.
Our air travel network simulations produced node-by-node seasonal SIS disease dynamics initiated by an outbreak of infection at a single airport. Examples are shown in Figure 2 for outbreaks at a highly connected node and a poorly connected node. A highly-connected Atlanta-seeded outbreak with β(t) = β 0 {(t) where {(t) is given by historical Google Flu Trends (see Methods and Figure 1(c)) shows a faster and higher amplitude dispersion of the disease compared to a poorly connected Colorado Springs-based outbreak. We chose a 50% initial infected fraction to magnify the properties of the network. The faster growth rate for Atlanta seeding was due to the number of flight connections and the higher peak amplitude of infection is due to the larger population of the Atlanta airport catchment area. The Atlanta seeded infection quickly moved to the other hubs in the network (the top seven hubs have larger dots) and remained concentrated in the South.
Due to the reduction of infectivity at seasonal minima and subsequent stochastic fadeout, the outbreak smoothly faded away for both cases within three years for this combination of β = 12 and γ = 3, or R 0 = 4. Again, this model was meant only to provide an example of an SIS outbreak rather than attempting to reproduce the fine structure of influenza dynamics. The expected contrast for a Colorado Springs-seeded outbreak is seen in the righthand panels (b,d). This outbreak at a less-connected node entered the network at the hubs and traveled more uniformly throughout the country, with a much lower intensity (note difference in scales for I(t)). These results indicate that our simulation was functioning as expected. Previously published simulations already noted that the connectivity of a node is significantly more important for spatial disease propagation than physical distance [21]. Short movies showing the evolution of these outbreaks are provided as electronic supplementary material.

Fractional reaction-diffusion with an outbreak on a line
The power law tail for the probability of displacements on the air travel network (Figure 1(a)) can be studied naturally with the non-local spacefractional diffusion operator (see Methods). In pursuit of a minimal model for non-diffusive transport on a realistic flight traffic network, we restricted the numerical solution of the fractional diffusion equation to one dimension. This avoided the problem of defining the fractional-derivative in two-dimensions and, as shown below, was sufficient for a quantitative agreement, in α and D α , for SIS disease dynamics between flight traffic and fractional-diffusion. This agreement was surprising since the one-dimensional fractional diffusion is structured very differently from the air travel network.
We implemented the SIS reaction system into this space-fractional reactiondiffusion model using both deterministic and stochastic solvers for the transport. As demonstrated previously [24], an operator splitting method with a multinomial kernel can efficiently and accurately solve the space-fractional diffusion equation for a range of α exponent values. These results are shown in Figure 3 for constant β = 12 and γ = 3. We chose two reference values of α only to highlight the differences between a very superdiffusive α = 0.3 and a less superdiffusive α = 1.3. The scaling exponent η of dispersion with time, σ 2 = D e t η , is more superdiffusive for smaller values of α: η = 2/α. The smaller value of α = 0.3 clearly leads to faster spreading of the disease compared to α = 1.3.
For the deterministic case, Figure 3(a,c), we used the solution of the fractional diffusion operator in Equations 6 and 7. By contrast, we introduced stochastic effects into the fractional-reaction-diffusion system using the Gillespie algorithm as described in the Methods. We used 482 nodes along a line as shown in Figure 3(c-d), seeded at node i = 241, with 1000 population units per node. The outbreak was seeded with all population units infected at the center node in order to magnify the amplitude and shorten the timescale for analysis. This choice only affects the dynamics quantitatively.
Stochastic behavior is apparent in Figure 3(b,d) that show the temporal and spatial dependence of the fractional SIS model for the same values of α in panels (a) and (c). We saw that stochasticity can affect the time of saturation for I i for nodes far from the center of the outbreak, though the steady-state mean fraction of infected is the same. The spatial dependence of I i on the grid position also shows the importance of stochastic reactions for amplifying infection at distal nodes, as well as the appearance of asymmetry due to fluctuations. Without stochastic fluctuations, the time delay to infection saturation is proportional to the distance from the outbreak node. Note that stochastic fluctuations have delayed the saturation time for the outbreak source node, especially in the less superdiffusive (α = 1.3/red) case.We chose to plot Figure 3c in log-log to highlight the shape of the propagating front, but chose lin-log for panel (d) to highlight the heterogeneity introduced by stochastic effects. Figure 3(a) shows the time-dependence of I i for three nodes taken from the symmetric solution shown in Figure 3(c). The number of infected population units decreased for the central outbreak node and rose for other nodes. One of the distal nodes is further from the outbreak, a difference that determined the delay in the arrival of the outbreak, separating the curves in Figure 3(a-b). All nodes eventually reached the same steady state with = 750 in the absence of seasonal forcing. Moving away from the center node delayed the arrival, as expected, of the outbreak wave. The shape of the propagating front of the infection in Figure 3(c) was determined partly by the fixed boundary condition, I bnd (t) = 0, but this is irrelevant due to the small fraction of population in the distal nodes as I(t) → I ∞ .
In addition to the obvious effect due to the distance from the outbreak node, the time to saturation was also influenced by the fractional diffusion parameters α and D α , such that larger (less superdiffusive) values of α imply a longer delay to full network saturation. Following this observation, we exploited the outbreak curve's dependence on the parameters of the fractional diffusion model to find the remarkable comparison described in the next section.

Capturing features of air travel outbreak with fractional diffusion
We explored the validity of using a fractional-diffusion model to approximate disease dynamics on the full flight network. For the most straightforward comparison, we excluded stochastic effects and seasonal forcing, since our main conclusions here are not sensitive to these choices. This caused the US air travel epidemic to saturate at I ∞ (β, γ), as we saw for the fractional diffusion model in the previous section. However, while the saturated state was determined by the SIS parameters, the evolution of I(t) was determined partly by the structure of the transport model, as we saw in both the previous sections. Thus, the shape of the network-averaged evolution curve should be a good measure of whether the simplified superdiffusion model can represent the full air travel network.
We compared the two SIS models using the cumulative difference, Θ (L2norm), between < I(t) > /N curves for 0 < t < t ∞ . Here, < I(t) > is the network averaged number of infected and t ∞ is the time at which the I ∞ value is reached. Figure 4(a) shows the direct comparison for a single R 0 (β = 12,γ = 3) over a range of fractional derivative α values with D α = 5. Better agreement is seen for smaller values of α, with α = 0.5 providing the best agreement. Figure 4(b) shows the consistency of the agreement between the two solutions for different sets of R 0 for the same values of α = 0.5 and D α = 5. The asymptotic value I ∞ = N (β−γ) β depends only on the SIS parameters. The shape of the curve is dependent on both R 0 and α, as Figure  4(a) shows.
In Figure 4(c-d), we show Θ for a parameter scan of α and D α (see Figure  1(a)). The Θ-surface defined by the parameter scan is similar for (c), where β = 12 and γ = 3 and (d) where β = 6 and γ = 3. Recall that R 0 = βS 0 γN . Regardless of R 0 , there is a unique minimum in α for each value of D α , consistent with the scaling in the long-tail limit of Equation 1. In other words, when α is larger, D α is forced by the form of q α to be larger with a similar scaling as observed in the Θ scan. While the value of D α is only an estimate, it is clear, at least, that fractional superdiffusion is qualitatively better than classical Brownian diffusion for this human mobility network.

Discussion
Devising a campaign strategy for the control and elimination of an infectious disease is a task that can benefit from computational modeling. We sought to encapsulate human mobility-driven spatial spreading of an infectious disease with a fractional reaction-diffusion model as an approximation to the full air travel network. We found good quantitative agreement considering that the mobility network can be characterized only partly by a power law tail. We conclude that a space-fractional diffusion model with superdiffusive statistics can emulate the dynamics of a seasonal disease, such as influenza, on a densely connected, long-range transportation network. This lightweight model offers another tool to explore the possible outcomes of disease outbreaks in silico.
Several caveats apply to the application of the fractional reaction-diffusion approach. In the limit of very few long-range connections, a heavy-tailed process is not a good approximation. Moreover, for a network of very few flights, fractional-reaction-diffusion is inaccurate and unnecessary. We are also aware that our translation of step sizes in miles between airports into a linear series of nodes does not have an obvious mapping function. However, it does seem to produce a compelling agreement. This work could be advanced by considering other modes of transportation such as private vehicles and foot traffic. It may also be useful eventually to consider a 2-D fractional-diffusion model, though the implementation of non-local effects in two dimensions is a challenging area of research [28,31].

Methods
In this work we employed two computational approaches to understand the problem of an oscillatory disease outbreak with spatial dispersion of infectious disease carriers, namely: 1) a stochastic fractional-diffusion-reaction equation and 2) a direct simulation of air traffic transport, the rates of which are derived from real-world flight networks. Both methods use the same algorithm for SIS infection kinetics separately at each node, including an option for stochastic kinetics. Our intent was to find empirically relevant parameter regimes where fractional-reaction-diffusion can reproduce macroscopic features of the spread of an infection on air traffic networks. All simulations were coded in Matlab and each time course (for 482 network nodes in the USA) required no more than 48 hours of serial calculations on a desktop computer. Computational time scaled linearly with the number of nodes as it is dominated by the SIS dynamics.

Air travel network simulation
We built a network reaction-diffusion model in Matlab that allows a great deal of flexibility. Transport can be specified with a matrix of transfer rates from node to node and reactions can be specified with a stochastic simulation algorithm, as described below. Source code is available upon request. For this paper, we included air travel in our simulations using flight schedules taken from the OAG airline database that includes flight data from over 900 airlines. This was the data source for all national flight networks used to build the transfer matrices. For each flight itinerary, the database includes origin and destination airport codes (IATA), number of seats and frequency of flight per week. We computed displacements of passengers summed over a week, assuming that all seats in the airplane are occupied. This approximation is sufficient for tracking infections for diseases with incubation times lasting several days.
Since we were interested in tractable simulations of stochasticity on the network of population centers served by air travel, we approximated each population with a proxy number of individuals, which we call population units. Each node contains this reduced number of units as a metapopulation, in contrast to a full agent-based model simulating each individual. Rather than trying to measure the relative size of each population served by an airport, we used the relative connectivity of each airport service center. Connectivity, C n , for a given node is the number of other nodes with a connecting flight to that node. A minimum node proxy population size is fixed at N min = 1000, and the total number of population units at a node, N tot = N min + B n C n , where B n scaled the connectivity to the total number of flights between the nodes. The population units in the simulation could be recalibrated for other applications.
As the population sizes were thus scaled down as described in the previous paragraph, the rate of transport from the flight network was renormalized to give relevant transport rates. In particular, the transfer rate of population was the sum i S i A i , where S i is the number of seats for flight i between the two nodes and A i is the number of days per week that the connection is serviced. Therefore, the transfer operator was applied on a weekly basis, which sets the time-scale for the simulation. The transfer matrix was also symmetrized to eliminate irregularities in the database in favor of a conservation of total network population. Moreover, the database-derived transfer operator was reduced by a numerical factor (F = 1 × 10 −7 ) to prevent the artificially small populations from being quickly depleted. Thus, our simulations on the air travel network should be considered as an approximation that preserved the relative size of each population center and the relative transport of population units between them. Applying the transfer operator on a weekly basis assumes that the infected proportion of the population does not change significantly on the weekly timescale.

Fractional diffusion simulation
A discrete-space, continuous time reaction-diffusion system can be modeled with the master equation, which is a differential equation [33]. The master equation describes the time evolution of a probability density vector, where each element of the vector denotes a state of the system. In population modeling, this state is typically the number of individuals in certain categories, such as infected with or susceptible to a disease. We are omitting some details of this method that can be found in other work [24], but including enough here for a degree of completeness.
We let n λ denote the number of people at a cell indexed by λ ∈ λ, n {n λ } λ∈λ , and p = p(n, t) the probability of finding n λ individuals at time t. Then, the master equation is [33]: where the operator A denotes diffusion, B denotes the transitions due to reactions. The master equation has an initial condition specified by p(n, 0) = p 0 (n), where p 0 (n) is typically a delta function. The formal solution to the master equation is: p(n, t) = e t(A+B) p(n, 0).
As discussed in [34] the master equation in some instances can be solved analytically, and therefore it is beneficial to partition the master equation and use analytical solutions when possible. To this end, we split the operators in the reaction-diffusion equation in order to use analytical solutions for the diffusive part, which yields [35]: which is equivalent to a simulation method where integration in time can be performed via the composition of the two operators: for independent samples s = 1, . . . , N. In other words, the diffusion process is simulated over a time-step ∆t and the system is updated, then the reactions are simulated with another time-step ∆t. These two steps constitute a single iteration in the numerical method. The reactions operator is discussed below. The diffusion operator φ consists of transitions that are not necessarily limited to adjacent cells. We define the fractional-diffusion kernel using a spatially second-order method (see Appendix, where the discretization is described in detail): where Γ(.) denotes the Γ function, D denotes the diffusion coefficient, G (λ ′ ) j denotes the jth component of the vector G (λ ′ ) , h the size of the compartment, α the value of the fractional derivative, and j G j (t) = 1 provided that the following time-step criterion is satisfied: In practice, the timestep should be chosen as an order of magnitude less than the right-hand side of Equation 8. The spatial boundary condition fixes the value of the population at zero. If U = U(x, t) denotes the concentration of people at position x ∈ R at time t, then the expected value of the method described above converges to the symmetric, space-fractional diffusion equation [36,37,38,39,34]: where α ∈ (0, 2] and D α is the effective diffusion coefficient.

Reaction Dynamics
There has recently been significant work on efficiently simulating spatial processes using the stochastic simulation algorithm (SSA) or Gillespie's algorithm [40,41]. In these models, the domain is discretized into cells such that each cell is assumed to be homogeneous. The result of the discretization is a set of random variables that denote the number of discrete quantities in each cell [42]. Diffusion events are modeled as transitions to adjacent cells, where the rate incorporates the macroscopic diffusion coefficient and the distance between the cells. The transition rates are derived by using finite-volume methods and the spatial process has been shown to converge to classical diffusion in the limit of a large numbers [43,34]. Stochastic simulation algorithms can be used to simulate the reactions that occur within cells due to the homogeneity assumption. Recently there has been a marked interest in developing efficient algorithms for the simulation of inhomogeneous systems with classical diffusion processes [44,45,34,42,46,47].
For our simulations of fractional diffusion and air traffic, each node has a finite population subject to infectious disease kinetics. We use the standard susceptible-infected-recovered (SIR) and susceptible-infected-susceptible (SIS) models that we implemented numerically with either a deterministic or a stochastic simulation. Conceptually, the infection kinetics are applied to the nodes with a reaction operator φ (∆t) B . Since we are generally interested in studying stochastic effects, such as fade-out and bursting, we considered various stochastic simulation algorithms, such as Gillespie's algorithm [40,41], τ -leaping [48], or R-leaping [49] (see [50] for a review). For an overview of the master equation formulation of stochastic kinetics, see Van Kampen [33]. We successfully tested the code with the SIR model system, we focus here on the standard SIS system, where the infection contact rate potentially varies with time: where I is the number of infected, S is the number of susceptible, N is the total population and β and γ are parameters that determine the reproduction The β and γ parameters affect the rate at which the infection grows and the steady-state infected fraction. The variation in infection rate can be modeled with a simple sinusoidal function, or as we have implemented, with a time-course explicity fit to the Google Flu Trends (GFT) database [51]. The GFT data gives an estimate of the occurrence of influenza-like illness (ILI) according to internet searches. While the predictive power of GFT data has recently been challenged [52], the qualitative pattern of the GFT data is descriptive of influenza occurrence. In simulations of the USA flight network, we use the historical GFT time course for individual states, shown in Figure  1(c).
We scanned through γ and β to determine whether it is possible to use the fractional-reaction diffusion equation to minimally model the epidemic on the air travel network. We note that R 0 in pandemics such as the 1918 flu has been estimated to be in the range of range 1.4 to 2.8 [53]. Simulations with values of R 0 closer to unity exhibit greater stochastic effects [54]. In our results, a value of R 0 = 4 or R 0 = 2 is used as a reference scenario without having qualitative influence on the result.

Acknowledgments
The authors thank Bill and Melinda Gates for their active support of this work and their sponsorship through the Global Good Fund. Thanks also to Alexandre Bovet, Daniel Klein, Edward Wenger and Josh Proctor for stimulating conversations.

Author Contributions
BSB and PAE conceived the project. KBG analyzed the flight data and performed the numerical simulations and applications. BSB derived and tested the numerical method. BSB, KG, and PAE wrote the manuscript.

Competing interest
We declare no competing interests with respect to this work.

Figure 3
One dimensional spatial network SIS simulations with R 0 = 4 using fractional diffusion. One thousand population units are simulated at each node Blue denotes a very superdiffusive α = 0.3 and red denotes a less superdiffusive α = 1.3 in all panels.  Figure 4 (a) Time traces for a constant R 0 = 4, SIS-modeled outbreak. The USA flight network (solid curve) is compared to discretized fractional diffusion (dashed lines). As the fractional diffusion exponent α decreases the transport becomes more superdiffusive and the dashed line is closer to the USA network result. The best-fit value of this exponent, α = 0.5, is similar to the one measured in Figure 1(a). Smaller values of α cause the curve to overshoot the air travel result. (b) Comparison between fractional diffusion and USA flight network outbreaks for superdiffusive α = 0.5 for several values of R 0 using SIS dynamics at each node. (c-d) Result of parameter scan through α and D α displaying the absolute distance Θ between the growth curves for average number of infected. The two panels show the difference for constant values of (c) R 0 = 2 and (d) R 0 = 4. These scans show that there is an optimal value of α that nearly matches the air travel model.

Movie 1
Outbreak of a seasonal, stochastic SIS infection in the Atlanta International airport (ATL) catchment area. The colorbar represents the number of infected population units at each node of the air travel network. Population sizes are scaled relatively to the size of the catchment area. The infection spreads quickly to hubs in the network and then fades out as the seasonal infectivity falls and stochastic fluctuations eliminate infected units. The average number of infected for this simulation is shown in Figure 3(c).

Movie 2
Outbreak of a seasonal, stochastic SIS infection in the Colorado Springs airport (COS) catchment area. The colorbar represents the number of infected population units at each node of the air travel network. Population sizes are scaled relatively to the size of the catchment area. The infection spreads quickly to hubs in the network and then fades out as the seasonal infectivity falls and stochastic fluctuations eliminate infected units. The average number of infected for this simulation is shown in Figure 3